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ABSTRACT 

The  research  described  in  this  report  has  investigated  methods  for  pre¬ 
dicting  optimal  mean  square  estimation  error  in  nonlinear  estimation  problems 
associated  with  passive  acoustic  tracking.  The  objective  was  to  develop  per¬ 
formance  prediction  methods  that  are  computationally  efficient,  applicable 
to  realistic  passive  tracking  models,  and  accurate.  Previous  work  extended 
Cramer-Rao  methods  to  obtain  a  method  that  was  computationally  efficient  and 
applicable  to  a  large  class  of  realistic  mathematical  models.  The  present 
work  reported  here  applies  this  method  to  study  the  effect  of  an  uncertain, 
unstable  source  frequency  and  the  effect  of  a  broadband  source  signal  on  pas¬ 
sive  tracking  using  omnidirectional  and  directional  sonobuoys. 

In  addition,  the  report  describes  work  on  developing  performance  predic¬ 
tion  methods  that  are  more  accurate  than  Cramer-Rao  methods  in  low  signal-to- 
noise  ratio  cases.  One  method  is  based  on  rate  distortion  theory  and  shows 
great  promise.  This  method  requires  no  simulation,  it  is  analytically  and 
efficiently  computable  for  a  large  class  of  static  nonlinear  problems,  and  it 
is  better  than  the  Cramer-Rao  method  when  signal-to-noise  ratio  is  low.  The 
work  reported  here  also  investigated  a  numerical  method  of  performance  pre¬ 
diction  based  on  classical  ambiguity  analysis.  This  investigation  obtained 
bounds  on  the  error  between  the  numerical  prediction  and  the  exact  mean  square 
error  as  a  function  of  the  size  of  the  numerical  computation. 
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SECTION  1 
INTRODUCTION 


1 . 1  INTRODUCTION 

The  operational  deployment  and  mission  objective  of  airborne  acoustic  sur¬ 
veillance  (passive  sonobuoys  monitored  by  aircraft)  create  a  special  acoustic 
environment  and  restrict  the  information  processing  capabilities  in  respects 
which  differ  from  other  passive  acoustic  surveillance  methods  (towed  arrays, 
hull-mounted  arrays,  bottom-anchored  arrays).  But  the  currently  operational 
and  planned  next  generation  of  airborne  signal  processing  systems  have  an 
architecture  similar  to  other  passive  acoustic  systems.  That  is,  the  proces¬ 
sing  system  consists  of  (1)  a  "front-end"  which  performs  spectral  analysis  of 
the  raw  signal  from  each  sonobuoy  to  extract  frequency  and  bearing  informa¬ 
tion  about  the  target  and  (2)  a  "back-end"  which  uses  the  front-end  output  to 
detect,  locate,  or  track  the  target.  Much  work  has  been  done  on  improving  the 
design  of  the  individual  modules,  but  there  appears  to  have  been  little  work 
on  evaluating  and  designing  the  overall  signal  processing-tracking  system. 

The  research  described  in  this  report  addresses  the  problem  of  evaluating 
the  potential  performance  that  might  be  achieved  by  improving  the  design  of 
the  overall  signal  processing  and  tracking  system.  Specifically,  our  research 
problem  is  to  develop  mathematical  methods  and  numerical  algorithms  to  estimate 
the  optimal  tracking  performance  possible  with  given  mathematical-physical 
models  of  acoustic  signals  and  sensors.  Our  objective  has  been  to  develop 
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performance  prediction  methods  that  are  computationally  efficient,  applicable 
to  realistic  passive  tracking  models,  and  accurate.  In  our  previous  work  [1]* 
we  developed  a  Cramer-Rao  method  to  obtain  a  method  that  was  computationally 
efficient  and  applicable  to  a  large  class  of  mathematical  models.  In  Section 
2  of  this  report  we  have  shown  that  this  method  is  easy  to  apply  to  more 
realistic  models  than  the  ones  used  in  [1].  Specifically,  we  have  used  the 
method  to  study  the  effect  of  uncertain,  unstable  source  frequency  and  the 
effect  of  the  presence  of  a  broadband  source  component  on  tracking  accuracy. 

In  some  nonlinear  estimation  problems  of  low  signal-to-noise  ratio, 
Cramer-Rao  methods  may  predict  performance  much  better  than  the  optimal  pro¬ 
cessing  algorithm  can  actually  achieve.  This  disadvantJ.  f  Cramer-Rao  meth¬ 
ods  motivated  us  to  investigate  performance  prediction  methods  which  would  be 
more  accurate  when  the  signal-to-noise  ratio  was  low,  but  which  are  still 
efficiently  computable  for  a  large  class  of  realistic  models.  Sections  3  and 
4  focused  on  this  problem. 

Section  3  investigated  an  analytical  (i.e.,  not  requiring  simulation) 
method  based  on  rate  distortion  theory  [4],  This  method  shows  great  promise 
because  it  is  efficient  to  compute  for  a  large  class  of  nonlinear  problems 
and  it  is  better  than  the  Cramer-Rao  method  when  signal-to-noise  ratio  is  low. 
However,  the  method  requires  further  development  to  make  it  applicable  to 
realistic  dynamic  problems. 

Section  5  investigated  a  numerical  method  of  performance  prediction  often 
described  as  ambiguity  analysis  [7], [15].  This  method  is  essentially  based 
on  numerical  computations  rather  than  on  analytical  formulas.  The  method  can 

*References  are  indicated  by  numbers  in  square  brackets,  the  list  appears  at 
the  end  of  the  main  body  of  this  report. 
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give  an  accurate  performance  prediction  provided  sufficient  computational 
resources  are  available.  Our  investigation  studied  the  relationship  between 
prediction  accuracy  and  computational  complexity  for  this  method.  Further 
work  remains  to  determine  the  precise  effect  of  signal-to-noise  ratio  on  the 
relationship  between  prediction  accuracy  and  computational  complexity. 

In  the  remainder  of  this  section  we  present  a  summary  of  the  research 
detailed  in  Sections  2,  3,  and  4. 

1.2  CRAMER-RAO  PERFORMANCE  ANALYSIS  OF  FREQUENCY  UNCERTAINTY 

AND  BROADBAND  SIGNALS 

In  our  previous  work  [1]  we  refined  available  Cramer-Rao  performance 
analysis  methods  to  exploit  special  features  of  the  airborne  acoustic  signal 
processing  and  tracking  problem.  In  particular,  it  was  possible  to  obtain 
an  efficient,  recursive  computation,  and  to  avoid  completely  Monte-Carlo  sim¬ 
ulation  despite  the  presence  of  nonlinear  measurements.  A  finite  dimensional 
stochastic  differential  equation  [2]  was  used  to  model  a  constant  velocity 
target  radiating  acoustic  signals  to  passive  omnidirectional  and  directional 
sonobuoys.  In  [1]  we  found  that  such  models  of  airborne  acoustic  tracking 
problems  have  the  following  structural  property:  all  nonlinearties  of  the 
state  equations  are  functions  only  of  the  source  kinematic  states  (position 
and  velocity).  The  other  state  variables,  modeling  such  quantities  as  source 
frequencies,  broadband  component,  received  signal  phase,  etc.  occur  linearly 
in  the  state  equations.  Because  of  this  structure  of  the  mathematical  model, 
we  were  able  to  derive  an  efficient  Cramer-Rao  type  lower  bound  that  treats 
the  source  kinematic  variables  as  unknown  parameters  and  all  other  state 
variables  as  random,  normally  distributed  parameters. 
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In  the  work  of  Section  2  we  have  illustrated  further  the  utility  of 
this  approach  by  analyzing  two  effects  not  studied  in  [1):  the  effect  of  an 
initially  uncertain  and  randomly  unstable  source  frequency,  and  the  effect 
of  a  broadband  component  of  the  source  frequency.  The  random  frequency  was 
parameterized  by  two  variables,  initial  root-mean-square  uncertainty  and  rate 
of  variation.  The  rate  of  variation  (studied  for  .01  to  1.0  Hz/min)  appeared 
to  have  negligible  effect  on  both  position  and  velocity  tracking  error. 

Initial  uncertainty  had  little  effect  on  position  tracking  error,  but  it  did 
have  a  significant  effect  on  velocity  tracking  error.  When  the  initial  source 
frequency  uncertainty  reaches  1  Hz,  the  initial  velocity  tracking  error  is 
not  substantially  reduced  until  the  source  passes  through  the  sonobuoy  field. 
This  indicates  that  initial  uncertainty  concerning  source  frequency  can  make 
the  velocity  tracking  performance  sensitive  to  source-sensor  geometry  (i.e., 
good  velocity  tracking  will  depend  more  crucially  on  good  geometry). 

The  broadband  source  component  was  modeled  as  a  simple  stationary  first- 
order  Markov  process.  The  bandwidth  of  this  process  was  fixed  at  200  Hz 
and  the  ratio  of  broadband  to  narrowband  power  was  varied  from  10“ 5  to  101*. 
Increasing  this  ratio  increases  the  total  source  power,  and  the  position  and 
velocity  tracking  error  decrease  as  a  consequence.  The  decrease  in  tracking 
error  is  comparable  to  the  decrease  in  error  with  increased  narrowband  source 
signal  power  studied  in  [1].  This  indicates  that  broadband  source  energy  is 
comparable  in  value  to  narrowband  source  energy,  and  that  the  broadband  com¬ 
ponent  of  the  source  signal  can  be  profitably  exploited  by  an  acoustic  signal 
processing  system. 

Many  other  realistic  models  can  be  analyzed  using  the  methods  of  Section 
2.  However,  before  analyzing  the  performance  of  other  realistic  models  of 
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passive  acoustic  tracking  problems,  we  need  to  determine  the  degree  of  opti¬ 
mism  inherent  in  the  performance  prediction  of  Section  2.  Sections  3  and  4 
present  one  approach  to  doing  this,  namely  by  trying  to  develop  more  accurate 
performance  predictions  with  which  to  compare  the  methods  of  Section  2.  It 
is  also  desirable  to  compare  these  performance  predictions  to  the  performance 
of  actual  algorithms.  The  methods  of  Section  2  suggest  a  processing  algorithm 
architecture  that  might  realize  the  performance  prediction  in  some  cases  (see 
[1]  for  discussion). 

1.3  RATE  DISTORTION  PERFORMANCE  ANALYSIS 

Communication  theory  provides  a  useful  interpretation  of  tracking  prob¬ 
lems  different  from  the  more  conventional  statistical  estimation  theory  point 
of  view.  Messages  are  generated  by  a  source  and  coded  by  an  encoder.  The 
encoded  messages  are  transmitted  through  a  channel,  decoded  by  a  decoder,  and 
received  by  a  user.  In  communication  problems  the  source,  channel,  and  user 
are  specified,  and  the  problem  is  to  design  encoder  and  decoder  so  that  mes¬ 
sages  received  by  the  user  are  accurate  reproductions  of  the  messages  gener¬ 
ated  by  the  source. 

One  can  interpret  a  tracking  system  as  a  type  of  communication  system  in 
the  following  way.  In  this  interpretation  the  message  generated  by  the  source 
is  a  set  of  target  parameters  (e.g.,  positions  and  velocities  at  a  given 
time).  The  encoder  for  passive  tracking  problems  does  nothing  to  code  the 
source  message.  In  active  tracking  we  can  control  the  encoder  to  some  extent 
(e.g.,  increase  signal  strength).  The  encoder  and  channel  for  the  tracking 
problem  represent  the  transformation  between  target  parameters  and  sensor 
outputs.  One  might  also  include  preprocessing  of  sensor  outputs  as  part  of 
the  channel  if  that  preprocessing  is  already  specified.  Finally,  the  decoder 
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for  the  tracking  problem  is  the  tracking  algorithm  which  provides  estimates 
of  target  parameters  to  a  user.  In  tracking  problems  the  source  (target 
model),  encoder  and  channel  (measurement  model),  and  user  are  specified,  and 
the  problem  is  to  design  a  decoder  (tracking  algorithm)  so  that  estimates 
received  by  the  user  are  accurate  reproductions  of  the  parameters  generated 
by  the  source. 

The  communication  theory  viewpoint  is  useful  because  it  allows  us  to 
apply  to  the  tracking  problem  techniques  of  information  theory  which  do  not 
exist  in  statistical  estimation  theory.  The  techniques  relevant  to  tracking 
performance  analysis  involve  rate  distortion  theory  [4]  first  developed  by 
Shannon  (5),  16].  Distortion  is  a  measure  of  the  average  error  between  the 
message  generated  by  the  source  and  the  decoded  message  received  by  the  user. 
In  a  tracking  problem  it  could  be  the  mean  square  error  in  the  tracking  algo¬ 
rithm's  target  parameter  estimation. 

Using  rate  distortion  techniques  and  some  simple  extensions  of  them  in 
Section  3,  we  showed  how  to  compute  analytically  rate  distortion  lower  bounds 
of  mean  square  error  for  static  nonlinear  estimation  problems  with  additive 
Gaussian  noise.  Specifically,  we  obtained  a  lower  bound  of  the  mean  square 
estimation  error  for  any  specified  component  of  the  state  vector.  We  showed 
that  the  rate  distortion  bound  is  asymptotically  tighter  than  the  Cramer-Rao 
bound  in  the  limit  of  low  signal-to-noise  ratio. 

Based  on  present  results,  the  rate  distortion  bound  offers  a  better 
approximation  of  mean  square  performance  in  the  high  measurement  noise  regime 
than  the  Cramer-Rao  bound.  Furthermore,  the  rate  distortion  bound  requires 
little,  if  any,  more  computation  than  the  Cramer-Rao  bound.  Thus,  the  rate 
distortion  bound  appears  to  complement  the  Cramer-Rao  method  in  the  nonlinear, 
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high  noise  regime  where  the  latter  bound  is  known  to  give  overly  optimistic 
approximations  of  the  true  mean  square  error.  However,  in  order  to  make  the 
rate  distortion  theory  useful  for  the  dynamic  nonlinear  estimation  problems 
of  tracking,  we  must  develop  our  current  results  in  two  significant  ways: 

1.  it  is  necessary  to  obtain  a  simple  rate  distortion  bound  in 
the  case  of  a  vector  state  and  a  general  vector  measurement; 

2.  it  is  necessary  to  derive  a  recursively  computable  bound  for 
dynamic  estimation  problems. 

Other  directions  for  further  investigation  exist  beside  these  two  neces¬ 
sary  extensions.  One  direction  would  extend  the  bounds  to  problems  for  which 
the  state  is  an  unknown,  non-random  parameter  (or  a  mixture  of  random  and  non- 
random  parameters).  In  [1]  we  found  that  a  large  class  of  tracking  problems 
can  be  modeled  by  a  state  process  which  consists  of  an  unknown  deterministic 
component  and  an  unknown,  Gaussian  distributed  random  component.  Rate  distor¬ 
tion  theory  for  nonstatistical  sources  (e  -  entropy  methods  [4])  may  allow  us 
to  derive  such  results. 

Another  direction  is  to  study  the  effect  of  architecture  constraints, 
such  as  preprocessing  of  measurements,  on  tracking  estimation  performance. 

We  investigated  this  problem  in  (1]  using  Cramer-Rao  methods.  A  rate  distor¬ 
tion  approach,  based  as  it  is  on  information  theory,  would  provide  a  more 
general,  more  accurate  method  of  analyzing  architectural  constraints. 

1.4  AMBIGUITY  PERFORMANCE  ANALYSIS 

Ambiguity  analysis  (17], [15]  Chapter  10)  is  an  attempt  to  understand  the 
global  nature  of  a  parameter  estimation  problem.  This  is  in  contrast  to  Cramer- 
Rao  methods  which  provide  a  more  local  analysis  of  estimation  performance. 

The  Cramer-Rao  lower  bound  on  mean  square  estimation  error  will  be  an  accurate 
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estimate  of  true  optimal  performance  provided  that  it  is  possible  to  acquire 
or  maintain  an  estimate  near  the  unknown  parameter  at  all.  The  local  problem 
(addressed  by  the  Cramer-Rao  method)  is  to  analyze  accuracy  given  acquisition; 
the  global  problem  (addressed  by  ambiguity  analysis)  is  to  analyze  the  acqui¬ 
sition  performance. 

The  ambiguity  method  approximates  the  mean  square  error  of  the  maximum 
likelihood  estimator  by  forming  a  weighted  sum  of  the  Cramer-Rao  lower  bound 
with  a  finite  number  of  discrete  errors.  The  weights  are  probabilities  asso¬ 
ciated  with  the  finite  hypothesis  testing  problem  of  choosing  one  of  a  finite 
number  of  regions  in  the  parameter  space.  The  regions  were  selected  so  that 
one  large  region  (proportional  to  in  size)  contained  the  true  parameter. 

The  rest  of  parameter  space  was  divided  into  smaller  regions  (proportional  to 
N-1  in  size).  We  showed  that  this  method  is  different  from  the  exact  mean 
square  error  by  a  terra  proportional  to  N-1.  Thus,  the  method  converges  to 
the  exact  mean  square  error  as  the  number  of  regions  increases,  and  the  error 
of  the  approximation  is  inversely  proportional  to  the  number  of  regions. 

Further  work  is  required  to  determine  how  the  magnitude  of  the  measure¬ 
ment  noise  (or  equivalently,  the  signal-to-noise  ratio)  enters  into  the 
approximation  error.  This  result  will  clarify  how  large  the  number  of  regions 
needs  to  be  for  a  given  signal-to-noise  ratio.  The  convergence  analysis  also 
needs  to  be  extended  to  the  general  case  of  vector  states  and  measurements 
and  to  the  case  where  a  measurement  process  is  observed.  The  order  of  the 
approximation  error  is  expected  to  remain  the  same  in  these  generalizations 
but  a  more  precise  idea  of  the  size  of  this  error  would  help  us  understand 
the  computational  feasibility  of  applying  the  ambiguity  analysis  method  to 
analyze  the  performance  of  complex  estimation  problems. 
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SECTION  2 

CRAMER-RAO  PERFORMANCE  ANALYSIS  OF  FREQUENCY 
INSTABILITY  AND  BROADBAND  SIGNALS 

2 . 1  INTRODUCTION 

In  our  previous  work  [ 1 ]  we  refined  available  Cramer-Rao  performance 
analysis  methods  to  exploit  special  features  of  the  airborne  acoustic  signal 
processing  and  tracking  problem.  In  particular,  it  was  possible  to  obtain  an 
efficient,  recursive  computation,  and  to  avoid  completely  Monte-Carlo  approx¬ 
imation  despite  the  presence  of  nonlinear  measurements.  In  this  section  we 
illustrate  further  the  utility  of  this  approach  by  analyzing  the  effects  of 
unstable  and  unknown  source  frequency  and  of  broadband  source  signals  on  the 
predicted  optimal  tracking  performance.  We  describe  the  mathematical  model 
we  have  used  in  subsection  2.2.  Subsection  2.3  presents  the  results  of  nu¬ 
merical  runs  and  subsection  2. A  presents  conclusions  based  on  these  results. 

2.2  MATHEMATICAL  MODEL  WITH  FREQUENCY  INSTABILITY  AND  BROADBAND  SOURCE 

A  finite  dimensional  stochastic  differential  equation  [2]  was  used  to 
model  a  constant  velocity  target  radiating  direct  path  acoustic  signals  to 
passive  omnidirectional  and  directional  sonobuoys.  This  model  has  the  gen¬ 
eral  form 

dx  =  f(x)dt  +  G  dw  (2-1) 


dy  =  h(x)dt  +  dv 


(2-2) 
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where  x  is  the  finite  dimensional  state  vector  and  y  is  the  finite  dimensional 
vector  of  measurements.  The  state  noise  and  measurement  noise  were  assumed 
to  be  independent  Gaussian  white  noise  processes.  Note  that  f  and  h  were  non¬ 
linear  functions  of  x  and  G  was  a  constant  matrix. 

The  state  vector  x  used  in  the  model  consisted  of  the  following  components. 


x  = 


-»■ 

V, 


Nb 


Nb 


(2-3) 


■>  + 

In  Eq.  2-3  the  expressions  x^,  x2  denote  the  two  orthogonal  components  of  the 
position  of  the  target  relative  to  some  fixed  position  (the  origin  at  0,0). 

The  expressions  v2  similarly  denote  the  components  of  the  target  velocity. 
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The  expression  ^  denotes  the  phase  of  the  narrowband  component  of  the  received 
acoustic  signal  at  buoy  number  k.  Similarly,  s^  denotes  the  broadband  compo¬ 
nent  of  the  received  acoustic  signal.  The  expression  f  denotes  the  transmitted 
source  frequency.  There  are  N  total  sonobuoys. 

The  measurement  vector  y  consisted  of  the  following  components. 


Nb 


(2-4) 


where  for  each  k  the  expression  y^  denotes  either 


yk  “  7k, om 


(2-5) 


if  sonobuoy  k  is  an  omnidirectional  buoy,  and 


rk,ora 
^k,d  1 
^k,d  2 


(2-6) 


if  sonobuoy  k  is  a  directional  buoy 
rectional  channel  signal  and  y 
signals . 


k,d  1 


In  any  case  yi^om  denotes  the  omnidi- 

y  denote  the  two  directional  channel 
k,d2 


The  target  model  assumed  constant  velocity  motion,  namely 


ALPHATECH,  INC 


dx^  =  v  ^  dt 
dx2  =  v2  dt 
dv1  =  0 
dv2  =  0 


(2-7) 


Note  that  Eq.  2-7  contains  no  driving  noise  on  the  right-hand  side.  This 
was  important  in  developing  an  efficient  method  to  compute  performance  as  de¬ 
scribed  in  [1],  The  initial  position  components  x  (0),  x^O)  and  the  initial 
velocity  components  v  (0),  v2(0)  were  treated  as  unknown  parameters  rather 
than  as  random  variables.  This  is  in  contrast  to  the  treatment  of  f,  4^  and 
sfc  as  random  parameters. 

The  acoustic  signal  radiated  by  the  target  was  assumed  to  consist  of 
two  parts:  a  narrowband  component  and  a  broadband  component.  The  narrowband 
component  was  modeled  as  follows.  The  source  frequency  f  satisfied  the  sto¬ 
chastic  differential  equation 


df  ■  -  af(f  -  f)dt  +  dwf  (2-8) 

where  Of  >  0,  f  and  the  variance  Of2  associated  with  wf  were  known  constants. 
The  initial  variance  of  f ( 0 )  was  assumed  to  be  the  steady  state  variance 
given  by 

- — -  Of  2 

(f(0)  -  TToT)2  - —  (2-9) 

2  Of 

where  the  overbars  (  )  denote  mathematical  expectations.  The  physical 

meaning  of  Eq.  2-8  is  that  the  source  frequency  f(t)  at  time  t  is  an  unknown 
random  variable  that  varies  randomly  about  the  constant  nominal  value  f 
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(f  =  100  Hz  in  our  numerical  examples).  Note  that  the  average  variation  of 
f  per  time  constant  is  of  order 

Of  •  (2-10) 

That  is,  Eq.  2-10  gives  the  average  speed  of  variation  of  the  source  frequency 
(it  has  units  Hz/min  in  our  numerical  examples). 

The  source  frequency  f  is  Doppler-shifted  and  drives  a  random  phase  equa¬ 
tion  given  by 


+  %,k  •  <2-*» 

In  Eq.  2-11  z  and  z  denote  the  coordinates  of  the  position  of  sonobuoy 
k,l  k, 2 

k,  and  c  is  the  assumed  constant  speed  of  sound.  The  noise  processes  w,^ 
were  assumed  to  have  the  same  variance,  o,j,2,  associated  with  all  of  them.  The 
model  also  allowed  the  possibility  of  correlation  between  w^^  and  w^j  for 
k  *  j.  A  constant  correlation  coefficient  p  (0  <  p  <  1)  was  assumed  for  all 
such  cases. 

The  physical  interpretation  of  Eq.  2-11  is  obtained  as  follows.  Consider 
a  signal  y(t)  with  phase  4>k(t).  specifically 

y(t)  -  sin  <tfc(t)  (2-12) 

If  one  assumes  that  the  phase  is  initially  unknown  (i.e.,  distributed  uniformly 
over  360  degrees),  then  y(t)  is  a  wide  sense  stationary  process  [ 3 J  with  total 
power  1/2  and  two-sided  power  spectral  density  at  frequency  u>  given  by 
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Ill  +  f  2  +  U)2 
4  D 


4  IT 


+  (f  +  to)2 
D 


(2-13) 


This  result  assumes  that  the  Doppler-shifted  frequency  is  fixed  or  varies 
slowly  relative  to  the  phase  variation.  The  expression  in  Eq.  2-13  represents 
a  spectrum  with  peak  at  frequency  f jj  and  width  of  order  o^2. 

The  broadband  component  of  the  received  signal  is  modeled  by  the  sto¬ 
chastic  differential  equation 


dsk  =  -  ^  •  sk  dt  +  dws>k 


(2-14) 


where  as  >  0  and  the  variance  os2  associated  with  wsk  are  known  constants. 

The  signal  sk(t)  is  a  stationary  process  with  two-sided  power  spectral  density 
at  frequency  <o  given  by 


1 

2* 


°s 


2 


(og2  +  a2) 


(2-15) 


total  power  o82/ 2as  and  bandwidth  of  order  cig. 

The  processes  $k  and  sk  are  combined  to  form  the  received  signal  process 
as  follows.  The  omnidirectional  channel  process  yk>om  satisfies  the  equation 

dyk,om  =  (sk  +  A  sin  4>k)dt  +  dvk>om  (2-16) 

where  the  constant  A  and  the  variance  oic>ora2  associated  with  vicj0m  are  assumed 
known.  Note  that  the  total  power  of  the  narrowband  component  of  yk,om  *s  A2/2 
compared  to  the  broadband  component  power  os2/2as  and  the  noise  component 
power  spectral  density  Ok,om2/2ff  (the  white  noise  dvk^om  has  infinite  total 


power  of  course). 
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1 


I  ft 


\ 


The  directional  channel  signals 
given  by  the  following. 


y  and  y  satisfy  similar  equations 
k,dl  k,d2 


d\,n  m  +  A  sln  • 


(;i  -  zk,i> 


dt  +  dv 


t1!  ‘  2k,l)2+  (S2'  Zk,2>: 


1/2 


k,dl 


dyk,d2 =  +  A  sin  • 


(*2  "  Zk,2) 


(*1  "  Zk,l)2+  ^2"  Zk,2^ 


(2-17) 


_  dt  +  dv 

1/2  k>d2 


(2-18) 


2 .3  NUMERICAL  EXAMPLES 

This  section  describes  the  numerical  examples  we  computed  to  study  the 
effects  of  unknown,  unstable  source  frequency  and  the  effects  of  a  broadband 
source  signal  on  predicted  optimal  tracking  performance.  Tracking  performance 
was  measured  by  the  root-mean-square  (rms)  errors  in  position  and  in  velocity 
versus  time  after  initial  contact.  The  target-sensor  geometry  assumed  in  each 
example  is  shown  in  Fig.  2-1;  nominal  test  parameters  are  shown  in  Table  2-1. 
Note  that  the  measurement  noise  for  omnidirectional  channels  was  chosen  so 
that  a  single  directional  buoy  within  5  kft  of  the  source  could  determine 
source  frequency  to  within  .2  Hz  and  bearing  to  within  4  degrees  based  on  1 
minute  worth  of  raw  data.  Figures  2-2  and  2-3  show  respectively  the  position 
and  velocity  tracking  performance  for  the  nominal  test  parameters.  Note  that 
the  rms  error  is  plotted  on  a  logarithmic  scale  in  these  figures. 
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Figure  2-1 .  Target  Sensor  Geometry 


TABLE  2-1.  NOMINAL  TEST  PARAMETER  VALUES 


TARGET 

Initial  x  position 

-15 

kft 

Initial  y  position 

0 

kft 

Initial  x  velocity 

1.2 

kft/min  (approx.  10  kt) 

Initial  y  velocity 

0 

kf t/min 

Initial  position  uncertainty 

30 

kft 

Initial  velocity  uncertainty 

1.5 

kft/min 

Total  time  of  contact 

25 

min 

SOURCE  SIGNAL 

Narrowband  center  frequency 

100 

Hz 

Initial  uncertainty  of  frequency 

.01 

Hz 

Frequency  variation  rate 

.01 

Hz/min 

Narrowband  component  line  width 

.1 

Hz 

Broadband  bandwidth 

200 

Hz 

Broadband  .'narrowband  power  ratio 

io- 

-5 

SENSOR 

Number  of  directional  buoys 

1 

Number  of  omnidirectional  buoys 

2 

Interbuoy  distance 

5 

kft 

Buoy  phase  correlation  coefficient 

0 

Omnidirectional  channel  noise 

.11 

Directional  channel  noise 

.03 
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Table  2-2  shows  the  cases  computed  to  study  the  effects  of  frequency 
uncertainty  and  frequency  instability.  As  shown  in  Table  2-2,  two  parameters 
measuring  initial  source  frequency  uncertainty  (given  by  Eq.  2-9)  and  average 
rate  of  variation  (given  by  Eq.  2-10)  were  varied  independently.  The  effect 
on  tracking  position  error  is  shown  in  Figs.  2-4  through  2-15;  the  effect  on 
tracking  velocity  error  is  shown  in  Figs.  2-16  through  2-27. 

Table  2-3  shows  the  cases  computed  to  study  the  effect  of  a  broadband 
source  component.  One  parameter,  the  ratio  of  total  broadband  to  total 
narrowband  power,  was  varied.  This  ratio  is  given  by 


The  effect  on  tracking  position  error  is  shown  in  Figs.  2-28  through  2-32; 
the  effect  on  tracking  velocity  error  is  shown  in  Figs.  2-33  through  2-37. 

2 .4  CONCLUSIONS 

Table  2-4  shows  the  effect  of  source  frequency  uncertainty  and  instabil¬ 
ity  on  the  minimum  tracking  position  error  achieved  during  contact  and  the 
velocity  error  given  at  the  same  time.  In  Figs.  2-4  through  2-15  this  error 
is  the  minimum  of  the  position  error  curves.  Note  that  the  position  error 
minimum  occurs  at  14  minutes  after  initial  contact;  the  closest  point  of 
approach  occurs  at  12.5  minutes  after  initial  contact.  As  Figs.  2-4  through 
2-15  and  Table  2-4  indicate,  the  initial  uncertainty  in  source  frequency  has 
a  small  effect  on  position  tracking  error;  the  rate  of  source  frequency  vari¬ 
ation  appears  to  have  virtually  no  effect.  Table  2-4  indicates  that  the  rate 
of  source  frequency  variation  has  no  effect  on  velocity  tracking  error;  but 
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the  initial  source  frequency  uncertainty  has  a  greater  effect  on  velocity 
error  than  it  has  on  position  error.  The  effect  of  initial  source  frequency 
uncertainty  on  velocity  error  is  shown  in  Figs.  2-16  through  2-27  where  the 
initial  uncertainty  varies  from  .01  Hz  to  10  Hz.  These  figures  show  that  the 
initial  velocity  uncertainty  (1.5  kft/min)  is  reduced  very  little  until  the 
source  passes  through  the  sonobuoy  field  when  initial  frequency  uncertainty 
is  high.  This  indicates  that  initial  frequency  uncertainty  could  make  the 
velocity  tracking  performance  sensitive  to  target-sensor  geometry  (i.e.,  good 
performance  will  depend  more  crucially  on  good  geometry). 

Table  2-5  shows  the  effect  of  a  broadband  source  signal  component  on  the 
minimum  tracking  position  error  achieved  during  contact  and  the  velocity  error 
given  at  the  same  time.  In  Figs.  2-28  through  2-32  this  error  is  the  minimum 
of  the  position  error  curves.  As  before,  the  minimum  occurs  at  14  minutes 
after  initial  contact.  Increasing  the  ratio  of  total  broadband  to  total  nar¬ 
rowband  power  decreases  position  and  velocity  error  as  indicated  in  Table 
2-5.  However,  this  effect  is  the  expected  consequence  of  the  total  increase 
in  signal  power  relative  to  the  background  noise  level.  That  is,  the  total 
narrowband  power  and  the  background  noise  levels  are  held  constant  in  the 
examples  described  here.  The  decrease  in  tracking  error  with  increased  broad¬ 
band  source  signal  power  is  comparable  to  the  decrease  in  tracking  error  with 
increased  narrowband  source  signal  power  studied  in  [1],  This  indicates  that 
broadband  source  energy  is  comparable  in  value  to  narrowband  source  energy, 
and  that  the  broadband  component  of  the  source  signal  might  be  profitably 
exploited  by  an  acoustic  signal  processing  system. 


ALPHATECH,  INC 


TABLE  2-2.  FREQUENCY  EFFECT  CASES  STUDIED 


Initial  Uncertainty  (Hz) 


Variation 

Rate 

(Hz/min) 


.01 

.1 

1 

10 

.01 

Fig.  2-4* 

Fig.  2-16** 

Fig.  2-5 

Fig.  2-17 

Fig.  2-6 

Fig.  2-18 

Fig.  2-7 

Fig.  2-19 

.1 

Fig.  2-8 

Fig.  2-9 

Fig.  2-10 

Fig.  2-11 

Fig.  2-20 

Fig.  2-21 

Fig.  2-22 

Fig.  2-23 

1 

Fig.  2-12 

Fig.  2-13 

Fig.  2-14 

Fig.  2-15 

Fig.  2-24 

Fig.  2-25 

Fig.  2-26 

Fig.  2-27 

TABLE  2-3.  BROADBAND  EFFECT  CASES  STUDIED 
Broadband  to  Narrowband  Power  Ratio 


10~5 

10"2 

1 

102 

104 

Fig.  2-28* 

Fig.  2-29 

Fig.  2-30 

Fig.  2-31 

Fig.  2-32 

Fig.  2-33** 

Fig.  2-34 

Fig.  2-35 

Fig.  2-36 

Fig.  2-37 

*  position  error  versus  time 

**  velocity  error  versus  time 
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TABLE  2-4.  FREQUENCY  EFFECT  ON  MINIMUM  TRACKING  ERROR 


Variation 

Rate 

(Hz/min) 


TABLE  2-5.  BROADBAND  EFFECT  ON  MINIMUM  TRACKING  ERROR 
Broadband  to  Narrowband  Power  Ratio 


10~5 

10~2 

1 

102 

lO4 

.40  kft* 

.40  kft 

.36  kft 

.22  kft 

.12  kft 

.08  kft/min** 

.08  kft/min 

.07  kft/min 

.04  kft/min 

.04  kft/min 

*  minimum  position  error 


Initial  Uncertainty  (Hz) 


.01 

.1 

1 

10 

.40  kft* 

.41  kft 

.42  kft 

.43  kft 

.01 

.08  kft/min** 

.11  kft/min 

.14  kft/min 

.15  kft/min 

.1 

.40  kft 

.41  kft 

.43  kft 

.43  kft 

.08  kft/min 

.10  kft/min 

.15  kft/min 

.15  kft/min 

1 

.40  kft 

.41  kft 

.43  kft 

.43  kft 

.08  kft/min 

.08  kft/min 

.15  kft/min 

.15  kft/min 

**  velocity  error  at  time  of  minimum  position  error 
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Figure  2-2.  Position  Error  for  Nominal  Case 
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Figure  2-3.  Vel 
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Figure  2-5.  Frequency  Effect  on  Position  Error 
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Figure  2-6.  Frequency  Effect  on  Position  Error 
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Figure  2-8.  Frequency  Effect  on  Position  Error 
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Figure  2-9.  Frequency  Effect  on  Position  Error 
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Figure  2-10.  Frequency  Effect  on  Position  Error 
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Figure  2-11.  Frequency  Effect  on  Position  Error 
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Figure  2-12.  Frequency  Effect  on  Position  Error 
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Figure  2-15.  Frequency  Effect  on  Position  Error 


ALPHATECH,  INC 


(uuu/3j)()  H0HM3 


2-27 


Figure  2-16.  Frequency  Effect  on  Velocity  Error 
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Figure  2-19.  Frequency  Effect  on  Velocity  Error 
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Figure  2-21.  Frequency  Effect  on  Velocity  Error 
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Figure  2-22.  Frequency  Effect  on  Velocity  Error 
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Figure  2-25.  Frequency  Effect  on  Velocity  Error 
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Figure  2-26.  Frequency  Effect  on  Velocity  Error 
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Figure  2-27.  Frequency  Effect  on  Velocity  Error 
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Figure  2-30.  Broadband  Effect  on  Position  Error 
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Figure  2-31.  Broadband  Effect  on  Position  Error 
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Figure  2-35.  Broadband  Effect  on  Velocity  Error 
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SECTION  3 

RATE  DISTORTION  PERFORMANCE  ANALYSIS 


3.1  INTRODUCTION 

This  section  discusses  the  rate  distortion  theory  approach  to  analyze 
mean  square  error  in  statistical  nonlinear  estimation  problems.  We  present 
here  preliminary  results  for  static  estimation  problems  and  compare  the  rate 
distortion  and  Cramer-Rao-Van  Trees  approaches.  Based  on  preliminary  results, 
the  rate  distortion  method  gives  a  computable  lower  bound  which  is  tighter 
than  the  Cramer-Rao-Van  Trees  lower  bound  in  the  regime  of  low  slgnal-to-noise 
ratio. 

This  section  is  organized  as  follows.  Subsection  3.2  presents  the 
necessary  background  in  communication  and  rate  distortion  theory.  It  also 
sketches  the  communication  system  approach  to  statistical  estimation  problems. 
Subsection  3.3  investigates  the  static  estimation  problem  with  a  scalar  state. 
Subsection  3.4  extends  this  to  the  vector  state  case. 

Finally,  subsection  3.5  concludes  the  section,  discussing  other  work  and 
directions  for  further  investigation. 

3.2  INFORMATION  THEORY  BACKGROUND 

3.2.1  Communication  Theory  Point  of  View 

Communication  theory  provides  a  useful  interpretation  of  tracking  prob¬ 
lems  different  from  the  more  conventional  statistical  estimation  theory  point 
of  view.  The  block  diagram  of  a  communication  system  is  shown  in  Fig.  3—1. 
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Messages  are  generated  by  a  source  and  coded  by  an  encoder.  The  encoded 
messages  are  transmitted  through  a  channel ,  decoded  by  a  decoder,  and  received 
by  a  user.  In  communication  problems  the  source,  channel,  and  user  are  spec¬ 
ified,  and  the  problem  is  to  design  encoder  and  decoder  so  that  messages 
received  by  the  user  are  accurate  reproductions  of  the  messages  generated  by 
the  source. 


Figure  3-1 .  Communication  System  Block  Diagram 

Figure  3-2  shows  how  one  can  interpret  a  tracking  system  as  a  type  of  commun¬ 
ication  system.  In  this  interpretation  the  message  generated  by  the  source  is 
a  set  of  target  parameters  (e.g.,  positions  and  velocities  at  a  given  time). 
The  encoder  for  passive  tracking  problems  does  nothing  to  code  the  source  mes¬ 
sage.  In  active  tracking  we  can  control  the  encoder  to  some  extent  (e.g.. 
Increase  signal  strength) .  The  encoder  and  channel  for  the  tracking  problem 
represent  the  transformation  between  target  parameters  and  sensor  outputs. 
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One  might  also  include  preprocessing  of  sensor  outputs  as  part  of  the  channel 
if  that  preprocessing  is  already  specified.  Finally,  the  decoder  for  the 
tracking  problem  is  the  tracking  algorithm  which  provides  estimates  of  target 
parameters  to  a  user.  In  tracking  problems  the  source  (target  model),  encoder 
and  channel  (measurement  model),  and  user  are  specified,  and  the  problem  is  to 
design  a  decoder  (tracking  algorithm)  so  that  estimates  received  by  the  user 
are  accurate  reproductions  of  the  parameters  generated  by  the  source. 


TARGET 

PARAMETERS 


TARGET 

PARAMETERS 


TARGET 

PARAMETERS 

ESTIMATES  R-1989 


Figure  3-2.  Communication  System  Interpretation  of  Tracking  System 
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The  communication  theory  viewpoint  is  useful  because  it  allows  us  to 
apply  to  the  tracking  problem  techniques  of  information  theory  which  do  not 
exist  in  statistical  estimation  theory.  The  techniques  relevant  to  tracking 
performance  analysis  involve  rate  distortion  theory  [4]  first  developed  by 
Shannon  [5], [6].  Distortion  is  a  measure  of  the  average  error  between  the 
message  generated  by  the  source  and  the  decoded  message  received  by  the  user. 
In  a  tracking  problem  it  could  be  the  mean  square  error  in  the  tracking  algo¬ 
rithm's  target  parameter  estimates.  An  important  question  about  a  communica¬ 
tion  problem  is:  given  a  source,  channel,  and  user,  under  what  conditions  is 
it  possible  to  design  an  encoder  and  decoder  that  reproduce  the  source  output 
for  the  user  with  an  average  distortion  that  does  not  exceed  some  specified 
upper  limit  D?  This  question  is  the  analog  of  the  tracking  performance  analy¬ 
sis  problem  we  are  interested  in:  given  target  and  measurement  model,  under 
what  conditions  is  it  possible  to  design  a  tracking  system  that  estimates 
target  parameters  with  an  average  error  that  does  not  exceed  some  specified 
upper  limit  D? 

Rate  distortion  theory  Is  able  to  answer  the  communication  system  ques¬ 
tion  in  a  precise  and  relatively  simple  way.  Associated  with  the  source  and 
user  is  a  function  R(D)  of  D,  called  the  rate  distortion  function.  Associated 
with  the  channel  is  a  quantity  C  called  its  capacity.  One  can  achieve  average 
distortion  D  if  and  only  if  channel  capacity  exceeds  R(D).  A  typical  rate 
distortion  function  is  shown  in  Fig.  3-3.  To  apply  rate  distortion  theory  to 
the  tracking  problem  we  need  to  find  the  rate  distortion  function  R(D)  asso¬ 
ciated  with  the  target  model  and  the  user-defined  error  criterion,  and  find 
the  capacity  C  associated  with  the  measurement  model.  These  quantities  will 
tell  us  that  an  average  tracking  error  smaller  than  D  can  be  achieved  only  if 
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R(D)  <  C.  Note  that  the  converse  statement,  that  R(D)  <  C  implies  we  can 
achieve  average  error  of  D,  is  not  true  in  tracking  because  we  cannot  control 
the  encoder  as  one  does  in  conventional  communication  systems.*  Thus,  the 
inequality  R(D)  <  C  imposes  a  lower  bound  on  the  average  distortion,  lower 
than  which  we  cannot  achieve  with  the  given  target  model  and  measurement 
model.  However,  it  may  be  possible  that  no  tracking  algorithm  is  able  to 
achieve  this  lower  limit.  Thus,  rate  distortion  theory  will  give  us  lower 
bounds  on  tracking  error,  and  we  will  need  to  study  the  tightness  of  these 
lower  bounds  as  a  separate  issue. 

Rate  distortion  theory  provides  techniques  for  computing  or  approximating 
R(D)  and  C  for  general  classes  of  sources,  users  (i.e.,  fidelity  criteria), 
and  channels.  We  will  describe  some  of  the  fundamental  results  of  rate 
distortion  theory  in  the  next  subsection  and  apply  these  results  to  tracking 
performance  analysis  in  subsequent  subsections. 

3.2.2  Rate  Distortion  Theory  Fundamentals  and  Estimation  Problems 

We  are  interested  in  tracking  problems  which  can  be  formulated  as  the 
following  type  of  statistical  estimation  problem. 

x(t+l)  =  Ax( t)  +  w(t)  (3-1) 

y(t)  *  h(x(t))  +  v(t)  (3-2) 

In  Eqs.  3-1  and  3-2,  the  variable  x(t)  is  the  state  at  time  t  we  desire  to 
estimate  given  measurements  up  to  that  time.  Note  that  the  state  evolves 
linearly  (Eq.  3-1),  and  we  assume  w(t)  and  v(t)  are  zero  mean  Gaussian  random 
vectors.  The  performance  analysis  problem  is  to  approximate  the  mean  square 


*For  very  special  systems  (e.g.,  linear  Gaussian  ones),  this  is  true 
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error  of  an  optimal  estimator.  In  previous  work  [lj  we  found  that  many  track¬ 
ing  problems  can  be  formulated  as  above,  and  this  particular  mathematical 
structure  simplifies  the  computation  of  the  Cramer-Rao-Van  Trees  performance 
bound  [1] »[7]  ,[8]  ,[9]  .  We  wish  to  see  here  whether  similar  simplifications 
occur  for  rate  distortion  theory  performance  bounds. 

Before  attempting  to  tackle  the  dynamic  problem  formulated  in  Eqs.  3-1 
and  3-2  we  will  study  the  static  problem 

y  =  h(x)  +  v  (3-3) 

where  x,  v  are  assumed  to  be  Gaussian  random  vectors  and  h  is  a  nonlinear 
function  of  x.  Our  approach  is  to  understand  the  general  static  problem  of 
Eq.  3-3  first.  We  can  then  write  the  dynamic  problem  of  Eqs.  3-1  and  3-2  as 
a  large  static  problem  and  try  to  exploit  the  recursive  structure  of  this  spe¬ 
cial  type  of  static  problem  to  obtain  an  efficient,  recursive  approximation 
of  the  minimum  mean  square  estimation  error.  In  this  report  we  will  consider 
only  static  problems;  we  will  discuss  dynamic  problems  in  a  subsequent  report. 

The  rate  distortion  function  of  a  memoryless  scalar  Gaussian  source  of 
mean  x  and  variance  Q  with  respect  to  the  squared-error  criterion  is  [4,  p.99] 


1 

R(D)  =  -  log 
2 


(3-4) 


The  capacity  C  of  a  channel  defined  by 

y  -  h(x)  +  v  (3-5) 

where  v  is  zero  mean  Gaussian  with  covariance  matrix  R  and  dimension  m  is 
given  by  the  mutual  information  I(y;x)  between  y  and  x. 
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C  -  I(y;x)  (3-6) 

We  cannot  compute  I(y;x)  for  general  nonlinear  h,  but  we  can  approximate  it 
as  follows. 

I(y;x)  -  H(y)  -  H(y|x)  (3-7) 

Equation  3-7  gives  the  mutual  information  in  terms  of  the  differential  entropy 
H(y)  and  the  conditional  differential  entropy  H(y|x).  Now  we  can  compute 
H(y|x)  exactly: 

m  i. 

H(y|x)  =  H(v)  =  -  log(2ire  [det  R]m)  (3-8) 

We  cannot  compute  H(y)  in  general,  but  we  can  bound  it  as  follows  [4], 

m  _L 

H(y)  <  -  log  (2 ire  [det  A]m)  (3-9) 

2 

In  Eq.  3-8  the  mxm  matrix  A  is  the  covariance  of  y: 

A  *  E( [y  —  E(y) ]  [y  -  E(y)]T)  (3-10) 

where  E(.)  denotes  mathematical  expectation  and  T  denotes  matrix  or  vector 
transposition. 

If  D  is  the  minimum  mean  square  error 

D  =  min  E((x  -  x(y))2)  (3-11) 

x 

where  the  mimimum  is  taken  over  all  estimators  x  based  on  the  measurement  y, 
then  rate  distortion  theory  [4]  tells  us  that 
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From  Eqs.  3-8  and  3-9  we  see  that 


1  detA 

c  ‘  5  l08(d^> 


Combining  Eqs.  3-4,  3-12  and  3-13  gives 


(3-13) 


1  Q  1  detA 

2  108  (i>  ‘  2  lo«  W 


or  equivalently, 


QdetR 

D  >  "  :: 

detA 


Thus,  from  Eq.  3-11  we  see  that 


,  QdetR 

E((x  -  x(y))2)  >  - 

detA 


(3-14) 


(3-15) 


(3-16) 


for  any  estimate  of  x  based  only  on  y.  Note  that  in  this  problem  the  covari¬ 
ance  of  y  can  be  written 

A  =  T  +  R  (3-17) 

where 

r  -  E(fh(x)  -  E(h(x))]  [h(x)  -  E(h(x))]T)  .  (3-18) 


Thus,  we  have 


E(x  -  x(y))2) 


QdetR 

>  - 

det  (T  +  R) 


(3-19) 


This  is  the  basic  result  of  rate  distortion  theory  we  will  use  in  the  follow¬ 
ing  subsections  to  analyze  the  minimum  mean  square  error  in  static  estimation 
problems . 
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3 .3  SCALAR  STATE 

3.3.1  Computation  of  the  Rate  Distortion  Bound 

In  this  subsection  we  consider  the  problem  of  estimating  a  scalar 
Gaussian  state  x  with  mean  x  and  variance  Q  given  the  vector  measurement 

y  =  h(x)  +  v  (3-20) 


where  v  is  a  Gaussian  random  vector  with  dimension  m,  0-mean,  and  covariance 
R.  In  subsection  3.2  we  found  the  basic  rate  distortion  bound  (RDB)  of  mean 
square  error: 


E((x  -  x(y))2) 


Q  •  detR 
>  det  (T  +  R) 


(3-21) 


where 

T  -  E([h(x)  -  E(h(x))]  [h(x)  -  E(h(x))]T)  (3-22) 

To  compute  the  RDB  requires  computing  T,  or  equivalently,  computing 

E(h(x))  (3-23) 

and 

E(h(x)h(x)T)  .  (3-24) 


The  expectations  Eqs.  3-23  and  3-24  are  taken  with  respect  to  a  Gaussian  dis¬ 
tribution  and  can  be  computed  in  closed  form  for  a  large  class  of  nonlinear  h. 
Note  that  we  utilized  this  fact  in  earlier  work  [1]  to  compute  Cramer-Rao- 
Van  Trees  bounds.  Indeed,  if  h(x)  is  a  sum  of  products  of  polynomials,  expo¬ 
nential,  and  sine  or  cosine  functions,  then  we  can  compute  the  expectations 
in  Eqs.  3-23  and  Eq.  3-24  in  closed  form.  This  is  true  also  if  x  is  a  vector. 
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The  computation  of  such  expectations  derives  from  the  basic  formula 


where 


E(  $(u,x) )  =  'i'(u) 

$(u,x)  =  exp(uTx) 

¥(u)  =  exp(u^5T  +  i  u^Qu) 
2 


(3-25) 


(3-26) 


where  u  =  (uj,  U2 ,  ...un)^  is  a  vector  of  complex  numbers,  and  x  is  a  Gaussian 
random  vector  of  dimension  n,  mean  x,  and  covariance  Q.  Consider  the  scalar 
case  n=l  for  example.  If  n  is  a  real  number  we  find 

-  1  o 

E(eUX)  =  eux  +  ”  u  Q  •  (3-27) 

If  u  "  i  (imaginary  number  /-T) ,  we  obtain 


E(e 


(3-28) 


which  gives  the  two  results 

-1 

-  Q 

E(sin  x)  =  sin  x  •  e  2  (3-29) 

-1 

-  ; Q 

E(cos  x)  ■  cos  x  e  2  ,  (3-30) 


Taking  derivatives  of  Eq.  3-26  with  respect  to  u  gives  us  an  expression  for 
E(xn) : 


E(xn) 


du 


-  1  „ 
eux  +  -  u2q 

2 


u»0 


(3-31) 


We  can  obtain  other  expectations  by  combining  these  operations 
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In  general,  we  consider  the  class  An  of  all  functions  f(x)  which  can  be 
written 


P 

f(x)  =  l  ck  •  x)  (3-32) 

k=l 


where  c^  are  complex  constants,  uk  are  constant  complex  vectors,  and  are 
differential  operators  of  the  general  form 


Dk 


3  k(l) 


3  k(2) 


8  k(n) 


For  example,  in  the  scalar  n=*l  case 


(3-33) 


xn  =  (J)n  •  <KO,x)  .  (3-34) 

3n 

Note  that  if  fj  and  f£  are  in  An,  then  so  are  fj  •  f£  and  fj  +  f2»  The  class 
An  also  contains  all  constant  functions. 

If  f(x)  is  given  by  Eq.  3-32,  then 

E(f(x))  =1  ck  D^Ufe)  (3-35) 

k=l 

is  the  closed  form  expression  for  the  expectation. 

Thus,  we  see  that  if  each  component  function  of  h(x)  is  in  An,  then  each 
component  of  h(x)  h(x)T  is  also  in  An  and  we  can  compute  the  expectations  in 
Eqs.  3-23  and  3-24  in  closed  form.  Furthermore,  we  can  approximate  any  non¬ 
linear  function  h  using  functions  in  An5  and  this  gives  us  a  method  for 
approximating  the  expectations  in  Eqs.  3-23  and  3-24  for  general  h.  In  the 
next  subsection  we  will  use  this  approach  to  compute  the  rate  distortion  bound 


in  some  examples 
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3.3.2  Examples  of  the  Rate  Distortion  Bound 
3. 3. 2.1  Linear  Measurements 

Suppose  Chat  h(x)  *  h  •  x  is  linear.  Let  us  compute  the  RDB  for 

y  =  h  •  x  +  v  .  (3-36) 

The  bound  for  Eq.  3-36  will  be  the  same  as  for 

y  =  h  •  x  +  v  (3-37) 

where 

h  =  yflr1  •  h  (3-38) 


v  =■  v 


(3-39) 


The  second  version  will  simplify  computation  because  the  covariance  of  v  is 
just  the  mxm  identity  matrix  Im*  The  main  quantity  to  compute  is 


~  T 

det  (T+  Im)  ■  det  (h  h  •  Q  +  Im) 


(3-40) 


This  is  easy  if  one  notes  that  the  determinant  of  a  matrix  is  the  product  of 
its  eigenvalues.  The  matrix 

h  hT  •  Q  +  Im  (3-41) 

has  one  eigenvalue 

hT  h  ♦  Q  +  1  (3-42) 


and  m-1  eigenvalues  1.  Thus,  the  RDB  is 

Q 

RDB  -  - 

fi^ST  •  Q  +  1 


(3-43) 
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In  terms  of  h  and  R  this  is  just 

RDB  =  (hTR_1h  +  Q-1)-1  (3-44) 

which  is  in  fact  the  minimum  mean  square  error  for  the  linear  problem. 

3. 3. 2. 2  Scalar  Nonlinear  Measurements 

Let  us  assume  that  y  and  v  are  scalar  random  variables,  and  let  us  con¬ 
sider  examples  of  nonlinear  h(x).  For  simplicity,  we  will  assume  x  *  0  in 
this  set  of  examples. 

h(x)  =  xn 

The  RDB  in  this  case  is 


RDB  »  (Bn  •  Qn-lR-1  +  q— 1 )— 1 


where 

(2n) !  (n!>2 


2n  •  n!  2n [(n) ! ]2 
2 

for  even  n  =  2,  4,  6,...  and 

(2n) ! 

Bn - 

2n  ♦  n! 


(3-45) 


(3-46) 


(3-47) 


for  odd  n  »  1,  3,  5,...  .  Figure  3-4  shows  the  dependence  of  this  bound  on  R 
for  fixed  Q  *  1.0  and  n  ■  1,  2,  3. 

h(x)  *  sin  x  and  h(x)  =  x  -  x^/6 


The  bound  for  sin  x  is 


Distortion  Bounds  for  h(x] 
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ROB  -  (^  [1  -  e  Q  XR  x+  Q  x)  4 


U-4B; 


The  bound  for  the  third  order  expansion  of  sin  x,  namely  x  -  -  is 


RDB  =  (U  -  Q  +  —  •  Q2]R_1  +  Q-1)-1  . 

12 


(3-49) 


Figure  3-5  shows  these  bounds  together  with  that  for  h(x)  =  x  versus  R  for 
Q  -1.0. 


3. 3. 2. 3  Identically  Distributed  Conditionally  Independent  Scalar  Measurements 


Suppose  that  we  take  N  measurements 


y(t)  =  h(x)  +  v(t) 


(3-50 


such  that  v(t)  are  independent,  scalar  Gaussian  random  variables  of  variance 
R.  Let 


r  =  E([h(x)  -  E(h(x))]2)  . 


(3-51 


Then  we  can  compute  the  RDB  for  Eq.  3-50  in  terms  of  Q,  T,  R  and  N.  The 
static  vector  measurement  problem  equivalent  to  Eq.  3-50  has  the  RDB  given  by 


Q  •  det  (R  •  I[q) 
det  (Tn  +  R  •  IN) 


(3-52 


where  Ij q  is  the  NXN  identity  matrix  and  I*n  is  the  NXN  matrix  which  has  all  of 
its  elements  equal  to  T.  Reasoning  as  in  subsection  3. 3. 2.1  we  can  compute 
the  determinant  in  the  denominator  of  Eq.  3-52.  Note  that  is  a  rank  1 
matrix  with  only  one  nonzero  eigenvalue,  namely 


Tr  r«i  =  N  •  r 


(3-53 
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where  Tr  denotes  the  trace  of  a  matrix.  Consequently,  +  R  •  1^  has  one 
eigenvalue  equal  to  N  •  T  +  R  and  N  -  1  eigenvalues  equal  to  R.  Thus,  the 
determinant  is  (N  •  T  +  R)  •  RN-1  and  the  RDB  is 

RDB  =  ( [ N  •  T  Q-l]  R-l  +  .  (3-54) 


Note  that  as  N+“  the  RDB  is  asymptotically  equal  to 


RDB 


QR 

Nr 


(3-55) 


3. 3. 2. 4  Nonidentical,  Conditionally  Independent  Scalar  Measurements 
Suppose  that 

y(l)  =  x  +  v(l)  (3-56) 

y(2)  =  x2  +  v(2)  (3-57) 


where  v(l),  v(2)  are  independent  0  -  mean,  Gaussian  random  variables  with 
variance  R.  Assume  x  has  0  mean  and  variance  Q.  The  variance  of  the  vector 
measurement  y  =  (y(l),  yv2))T  is 


r 


Q  0 

0  2  Q2 


(3-58) 


Thus,  the  RDB  for  this  problem  is 

RDB  -  ([1  +  2Q  +  2Q2R-1]  •  R~1  +  Q"l)-1  (3-59) 

Figure  3-6  shows  this  bound  with  Q  =  1  versus  different  values  of  R.  Figure 
3-6  also  shows  the  mean  square  error  for  only  one  measurement  (namely  Eq.  3-56). 


Distortion  Bounds  for  One  and  Two  Measurements 
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3.3.3  Comparison  to  the  Cramer-Rao-Van  Trees  Bound 

3. 3. 3.1  The  Cramer-Rao-Van  Trees  Bound 

Before  comparing  the  Cramer-Rao-Van  Trees  lower  bound  (CRVB)  to  the  RDB, 
let  us  review  what  the  CRVB  is  for  the  static  estimation  problem  formulated 
above.  Van  Trees  [7J  has  derived  a  lower  bound  on  the  error  convariance  of 
any  estimator  x(y)  for  the  problem  in  Eq.  3-20.  The  bound  is 

E([x  -  x(y) ]  [x  -  x<y) ]T)  >  [T*  +  (3-60) 

where  the  inequality  is  in  terms  of  symmetric  matrices,  and  T*  is  defined  as 

3h  „  ,  3h 

T*  =  E( —  <x)Tr-1  — (x))  .  (3-61) 

3x  3x 

If  x  is  a  scalar  random  variable,  then  Eq.  3-60  gives  the  following  lower 
bound  on  the  error  variance 

E([x-i(y)]2)  >  IT*  +  .  (3-62) 

Note  that  Eq.  3-61  involves  computations  similar  to  those  required  for  the 

RDB.  Indeed,  if  h  is  a  member  of  the  class  of  functions  An  defined  in  subsec- 

3h 

tion  3.3.1,  then  the  components  of  —  also  belong  to  An.  In  this  subsection  we 

3x 

will  study  the  relation  of  the  CRVB  of  Eq.  3-62  to  the  RDB  of  Eq.  3-21.  Both 
are  lower  bounds  of  the  minimum  mean  square  error.  Can  we  determine  condi¬ 
tions  under  which  one  is  a  tighter  bound  than  the  other? 

3. 3. 3. 2  Comparison  of  Bounds  for  Scalar  Measurements 

Let  us  start  by  computing  the  CRVB's  corresponding  to  the  examples  of 


subsection  3. 3. 2. 2 
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h(x)  »  xn 

The  CRVB  in  this  case  is 


CRVB  =  (CnQn-l  R_1  +  Q-1)"1 


(3-63) 


where 


(2[n-l ] ) ! 
2n_1  [ n- 1 ] ! 


(3-64) 


Figures  3-7a  and  3-7b  show  CRVB  and  RDB  for  n  =  2,3  versus  R  with  Q  =  1.0. 
Note  that  ROB  >  CRVB  in  these  examples.  Indeed,  one  can  see  that 


Cn  >  Bn 


(3-65) 


for  all  n  =  1,2,...  and  therefore  CRVB  <  RDB  for  all  n. 


h(x)  «  sin  x  and  h(x)  =  x  -  x^/6 
The  CRVB  for  h(x)  ■  sin  x  is  found  to  be 

CRVB  =  (-  [1  +  e~2Q]R~l  +  Q”1)-1 
2 

x^ 

and  that  for  h(x)  »  x  -  ■ —  is 

6 


(3-66) 


3Q2  ,  .  . 

CRVB  =  ([1  -  Q  +  - ]R_1  +  Q-1)-1 

4 


(3-67) 


Figures  3-8a  and  3-8b  show  the  CRVB  and  RDB  for  these  two  nonlinear  examples 
(Q  *  1.0  and  R  is  varied).  Note  that  the  RDB  is  always  tighter  (i.e.,  CRVB  < 
RDB),  and  in  fact  one  can  prove  this  is  true.  It  is  interesting  to  note  that 
as  Eq.  3-48  predicts  (correctly)  that  the  mean  square  error  of  an  esti¬ 

mate  of  x  given  y  =  sin  x  +  v  blows  up.  The  CRVB  Eq.  3-66  predicts  (incor¬ 
rectly)  that  the  error  remains  bounded. 


Figure  3-7a.  Comparison  of  Rate  Distortion  and  Cramer-Rao-Van  Trees 
Bounds  for  h(x)  =  x2 


Figure  3-8b.  Comparison  of  Rate  Distortion  and  Cramer-Rao-Van  Trees 
Bounds  for  h(x)  ■  sin  x 
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General  nonlinear  h(x) 

As  one  might  expect  from  the  examples  above,  it  is  possible  to  show  that  CRVB 
<  RDB  all  the  time  in  case  of  scalar  measurements.  Furthermore,  one  can  show 
CRVB  =  RDB  if  and  only  if  h(x)  is  linear.  Recall  that  our  underlying  assump¬ 
tions  at  this  point  are  that  x  is  a  scalar  Gaussian  random  variable  and  y  is  a 
scalar  measurement.  We  prove  the  following  theorem: 

Theorem.  If  h(x)  is  continuously  differentiable  in  x,  and 
the  expectations  E(h(x)2),  E(Ih'(x)]2)  are  both  finite,  then 

CRVB  <  RDB 

where  CRVB  =  RDB  if  and  only  if 

h(x)  *  ax  +  b 


for  constants  a,b 

We  will  prove  ^his  theorem  in  the  remainder  of  the  subsection.  Assume 
first  that  x  =*  0  and  define 


<f(x) 


[h(x)  -  h(0)]2  x_1  e 


2Q 


(3-68) 


for  x  *  0,  and 


<S(0)  -  0 


(3-69) 


Note  that  *(x)  is  continuously  differentiable  and 

2 

x 

$'(x)  =  2lh(x)  -  h(0)]  h'(x)x-1  e"  2Q 
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for  x  *  0,  and 


r(0)  =  [h"(0)]2  . 


Given  that  E(h(x)2)  and  E(h'(x)2)  are  finite,  we  have  that 


/  4>'(x)  dx  *  lim  [<t(x)  -  $(-x)]  =  0 

—  00 

x+» 


Note  that  by  definition  of  mathematical  expectation. 


E(f (x) ) 


V211Q 


/  f(x)  e  2Q  dx  . 


Using  Eqs.  3-70  and  3-73  to  rewrite  Eq.  3-72  gives  us 


E(h'(x)2)  =  Q-l  E([h(x)  -  h(0)]2) 
+  E  ( [h'(x)  - 


Let  h  =  E(h(x)).  Then 


E(h'(x)2)  -  Q-l  E([h(x)  -  h] 2)  +  q-1  [h  -  h(0)]2 

h(x)  -  h(0) 


+  E([h'(x)  - 


)2) 


it  follows  that 


E(h'(x)2)  >  q-1  E([h(x)  -  h]2) 


with  equality  if  and  only  if 
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The  Eqs.  3-77  and  3-78  are  equivalent  to  h  being  of  the  form 

h(x)  =  ax  +  b  (3-79) 


for  constants  a,  b.  Note  that  Eq.  3-76  is  equivalent  to 


r**R  >  r*Q~l 


Since 


CRVB  “  (r*  + 


and 


RDB  *  (rQ"lR_1  +  Q-1)-1 


(3-80) 


(3-81) 


Eq.  3-80  proves  that  CRVB  >  RDB,  at  least  for  the  case  x  =  0. 

The  general  case  of  x  H  follows  easily  from  the  x  =*  0  case.  Simply 
apply  the  earlier  results  to  x  -x  with  the  measurement  function  h(x  +  x) . 
This  problem  will  yield  the  same  bounds  as  for  x,  h(x). 


3. 3. 3. 3  Comparison  of  Bounds  for  Identically  Distributed  -  Conditionally 
Independent  Measurements 

Under  the  problem  assumptions  of  subsection  3. 3. 2. 3,  we  can  show  that 
CRVB  <  RDB  holds  for  general  nonlinear  h(x).  Recall  (Eq.  3-54)  that 

RDB  =■  (  [N  *rQ-l  JR'1  +  q-l)“l  (3-82) 


One  can  easily  show  that 

CRB  -  ([N*r*R]R~i  +  Q-I)-I  (3-83) 

We  proved  that  rQ“l  <  T*R  in  the  last  subsection.  Consequently,  we  see  that 
CRVB  <  RDB  in  this  case  also.  Note  that  CRVB  predicts  a  mean  square  error 
that  is  asymptotic  to 


3-28 
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i 

CRVB  -  -  (3-84) 

Nr* 

as  N+®.  This  is  in  constrast  to 

QR 

RDB  *  —  .  (3-85) 

Nr 

Thus,  we  have 

CRVB  r 

-  -  -  (3-86) 

RDB  QRT* 

asymptotically  as  N+«.  For  nonlinear  h(x)  we  saw  previously  that  r/QRT*  <1. 
Thus,  the  CRVB  predicts  a  faster  rate  of  decrease  in  mean  square  error  than  is 
in  fact  possible. 

3. 3. 3. 4  Comparison  of  Bounds  for  Vector  Measurements 

From  the  preceding  results  one  might  conjecture  that  CRVB  <  RDB  in 
general.  The  following  example  shows  that  this  need  not  be  true.  Consider 
the  example  of  subsection  3. 3. 2. 4.  The  CRVB  for  this  example  is 

CRVB  -  ((1  +  4Q]R"1  +  Q“l)-1  (3-87) 

Figure  3-9  shows  CRVB  and  the  corresponding  RDB  of  Eq.  3-59,  namely 

RDB  =  ([1  +  2Q  +  2q2r-1]R"1  +  Q“l)~l  .  (3-88) 

It  is  easy  to  see  in  this  example  that  CRVB  >  RDB  if  R  <  Q  and  CRVB  <  RDB  if 
R  >  Q.  We  can  show  in  general  that  CRVB  <  RDB  if  R  is  sufficiently  large 
(i.e.,  as  signal-to-noise  ratio  approaches  0). 

Without  loss  of  generality  (by  using  the  same  transformation  as  in 
subsection  3. 3. 2.1)  we  can  assume  that  the  measurement  noise  covariance  is  a 
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scalar  constant  R  multiplied  times  the  mxm  identity  matrix  Im.  Then  we  have 

RDB  =  Q  *(det  [  Injxm  +  R-1  r])-l  (3-89) 

and 


m 

CRVB  =  ([[  E(h'(x)2)  ]R-1  +  Q-1)-1  (3-90) 

k=l  k 


where  hk(x)  are  the  components  of  h(x).  The  determinant  in  Eq.  3-89  can  be 
expanded  in  powers  of  R~*  so  that 

det  (Ijnxm  +  R"1  0  =  1  +  TrlR_1r]  +  0(R~2)  (3-91) 

where  0(R”2)  denotes  terms  of  order  R~2  or  higher  powers  of  R~1 .  The  trace  is 

n 

Tr[R_1r]  =  R_1  l  E((hk(x)  -  E(hk(x))]2)  (3-92) 

k=l 

Thus,  we  see  that 


n 

RDB  =  (  l  E(lhk(x)  -  E(hk(x))]2)  •  Q-l  VT1  +  Q_1 
k=l 

+  0(R-2))”1  (3-93) 


The  theorem  of  subsection  3. 3. 3. 2  implies  that 

A 

E([hk(x)  -  E(hk(x))]2)  Q“  <  E(hk(x)2)  (3-94) 

for  each  k  with  equality  if  and  only  if  hk  is  linear.  Consequently,  'f  h  is 
not  linear,  then  CRVB  <  RDB  for  R  suffiently  large.  If  h  is  linear,  then  both 
bounds  are  equal  to  the  minimum  mean  square  error. 


3-31 
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3 .4  VECTOR  STATE 

3.4.1  Rate  Distortion  Bound  for  Vector  State  and  Scalar  Measurement 

Rate  distortion  theory  naturally  gives  bounds  on  scalar  errors.  For 
vector  state  estimation  problems,  however,  we  need  a  bound  on  the  error  covar¬ 
iance  matrix  such  as  the  CRVB  provides  (Eq.  3-60).  In  this  subsection  we 
derive  a  RDB  for  vector  state  estimation  problems  with  scalar  measurements. 
Suppose  that 

y  -  h(x!,  X2)  +  v  (3-95) 

where  v  is  a  0-mean  Gaussian  random  variable  of  variance  R  and  x^ ,  X2  are 
jointly  Gaussian  random  variables.  We  are  interested  in  deriving  a  lower 
bound  for  the  mean  square  error 

E( fxi  -  xj(y)]2)  (3-96) 

of  an  estimate  xj(y)  of  xj  .  Suppose  X2  were  a  fixed,  known  constant.  Then 
the  previous  result  for  a  scalar  state  implies  that 

E([X1  -  xi (y) ]2 |x2)  >  ([r1(x2)Qr1]R“1  +  Qr1)’1  (3-97) 

where 

rl(x2)  *  E((h(xi,x2)  -  E(h(xj ,X2) |x2)]2|x2)  (3-98) 

and 

Q]  -  E( Ixj  -  E(xj  Jx2)]2|x2)  .  (3-99) 

Note  that  the  conditional  variance  in  Eq.  3-99  does  not  actually  depend  on  X2 
because  xj  and  X2  are  jointly  Gaussian. 


Equation  3-97  can  also  be  written  as 
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E((xi  -  xi(y)j2|x2)_1  <  i ri(x2)Qi~1jR_1  +  Ql"1 
Jensen's  inequality  states  that 

<J»(E(0)  <  E(<J>U)) 

if  f  is  a  convex  function  [10].  Note  that  $(  £)  =  (T1  is  convex  if 
Thus,  we  can  apply  Jensen's  inequality  to  obtain 

Eur1  <  Etc1) 

If  5  *  E([xj  -Xj(y)]2|x2)  we  obtain 

E([X1  -  x1(y)]2)~l  <  [Ti  •  Qi"1 ]R_1  +  Qi-1 

where 

Ti  =*  E(Ti(x2)) 

This  gives  the  rate  distortion  bound 

E( [x]  -  xi(y)]2)  ?  ( I TiQi- 1  JR-1  +  Ql-1)-1 

where 

Tj  *  E([h(xj,X2)  -  E(h(xj ,X2> |x2)]2) 


Note  also  that 

Ti  *  E(h(xi,x2>2)  -  E([E(h(xi ,X2) |x2>l2) 


Equation  3-105  is  our  basic  RDB  for  the  vector  state,  scalar 
case.  Note  that  X2  could  be  a  Gaussian  random  vector.  Thus,  for 

y  =  h(xi,X2,...,xn)  +  v 


(3-100) 

(3-101) 
£  >  0. 

(3-102) 

(3-103) 

(3-104) 

(3-105) 

(3-106) 

(3-107) 

measurement 

(3-108) 
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one  has  Che  same  bound  Gq.  3-105  except  Chat 

1*1  *  E(  [h(xj  . . xn)  -  E(h(xi,X2 . xn)  |x2...xn)]2)  .  (3-109) 

3.4.2  Computation  of  the  Rate  Distortion  Bound 

To  compute  the  RDB  of  the  previous  subsection  it  is  necessary  to  compute 
the  expression  Tj  in  Eq.  3-106  and  3-107  (or  more  generally,  Eq.  3-109).  It 
is  possible  to  do  this  in  much  the  same  way  as  we  did  in  subsection  3.3.1. 
Specifically,  if  h  belongs  to  the  class  An  of  functions  defined  in  subsection 
3.3.1  (Eq.  3-13),  then  we  can  compute  in  closed  form. 

Suppose  that  xj  is  a  random  n-dimensional  vector,  X2  is  a  random 
m-dimensional  and  xj ,X2  are  jointly  Gaussian.  Then  we  know 

E(xi ,X2)  -  a  +  Bx2  (3-110) 

E([xi-E(xi |x2)]  [xi  -  E(xi|x2)]T|x2)  -  C  (3-111) 

for  constant  vector  a  and  constant  matrices  B,  C.  If 

#(u,xj)  »  exp(uTxj)  (3-112) 

then 

E(*(u,xj)  jx2)  -  «u,x2)  (3-113) 

where 

T  T  IT 

iKutx_)  “  exp(u  a  +  u  Bx„  +  _  u  Cu)  .  (3-114) 

2 

where  u  is  an  n-dimensional  vector  of  complex  constants.  Note  that  for  any 
constants  a,  B,  C,  u  the  function  ip(u,x2)  x2  belongs  to  An.  Thus,  if 

P 

f(*l>  “  I  Ck  Djt#(ujc,x| ) 
k-1 


(3-115) 
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as  in  subsecdon  3.3.1  (Eqs.  3-32  and  3-34),  then 

P 

E(f (xj ) |x2)  -  l  Cfc  Dfc  lKuic^)  •  (3-116) 

k=  1 

In  other  words,  if  f(xj)  belongs  to  An,  then  E(f(xi)|x2)  belongs  to  A^. 

Consequently,  if  h(xj ,x2 , . . .xn)  belongs  to  the  class  An  of  functions, 
then  E(h(xj ,x2 . . . ,xn) |x2 , . . .xn)  belongs  to  An_i  and  so  does  [E(h(xi ,X2 , . . .xn) | 
X2,...xn)]2.  Thus  in  Eqs.  3-106,  3-107,  and  3-109  can  be  computed  in 
closed  form. 

3.4.3  Examples  of  the  Rate  Distortion  Bound 
3. 4. 3.1  Linear  Measurement 

Suppose  that  h(x)  =  h*x  is  linear  (h  is  a  row  vector  and  x  is  a  column 
vector).  Suppose  that  x  has  been  partitioned  into  a  one-dimensional  component 
X}  (which  we  want  to  estimate)  and  an  (n-1)  dimensional  component  x2. 

y  **  hj  *xj  +  h2  *X2  +  v  (3-117) 

is  the  measurement  equation.  The  mean  and  covariance  of  x  are  given  by 


E(x]  ) 

“  xi 

(3-118) 

E(x2) 

"  x2 

(3-119) 

E(  (xj 

-  xi 

)2)  - 

Qll 

(3-120) 

E(  (xj 

-  XI 

1  [x2 

-  x2]T) 

"  Ql2 

(3-121) 

E((x2 

-  x2 

]  Ixj 

-  xj ])  > 

"  Q21 

(3-122) 

E((x2 

-  x2 

1  (x2 

-  X2JT) 

-  q22 

(3-123) 
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The  condicional  mean  and  variance  of  xj  given  X2  are 

E(x1|x2)  -  X!  +  Q12  Q22-1  fx2  “  *2)  (3-124) 

E(l*l  -  E(xi  |x2)]21x2)  “  Qll  -  Ql2  Q22-1  Q21  •  (3-125) 

Thus,  we  find  that 

Ql  =  Qll  ”  Ql2  Q22-1  Q21  (3-126) 

and 

T  =  r2(x2)  =  hj2  •  Qj  .  (3-127) 

The  corresponding  bound  is 

RDB  -  (h^R-l  +  Qj"1)-1  (3-128) 

=  (hj2R-1  +  [Qjj  — Q12  Q22-1  Q211"1)'1  (3-129) 


Note  that  the  choice  of  the  component  x2  is  somewhat  arbitrary,  and  one  could 

try  to  select  it  to  make  RDB  as  large  as  possible.  For  example,  one  might 

T 

choose  x2  independent  of  xj  so  that  Q21  ■  Q12  ■  0.  Thus,  the  largest  RDB 
bound  obtained  by  choosing  x2  independent  of  xj  is  in  the  linear  case 

RDB  -  (hi2R-1  +  On-1)**1  .  (3-130) 

This  is  generally  smaller  than  the  minimum  mean  square  error.  We  will  examine 
this  more  closely  in  subsection  3.4.4  where  we  compare  the  RDB  with  the  CRVB. 
Let  us  remark  that  it  is  possible  to  develop  a  tighter  rate  distortion  bound 
that  gives  the  minimum  mean  square  error  exactly.  However,  this  bound  appears 
to  be  difficult  to  compute  in  nonlinear  problems. 
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3. 4. 3. 2  Nonlinear  Measurement 

Consider  the  simple  example 

y  *  xj2  +  *2  +  v  ,  (3-131) 

where  xj  and  X2  are  independent  Gaussian  random  variables  with  variances  Qn 
and  Q22  respectively.  Then  Qi  ■  Qll»  I*|  ■  2Qji2  and  the  RDB  is 

RDB  -  (UQulR-l  +  Qn"1)"1  .  (3-132) 

3.4.4  Comparison  to  Cramer-Rao  Van  Trees  Bound 
3. 4. 4.1  Examples 

The  following  examples  show  that  neither  CRVB  <  RDB  nor  RDB  <  CRVB  in 
general  for  vector  states  and  scalar  measurements.  Suppose  that  X}  and  X2  are 
Independent  Gaussian  random  variables  with  respective  variances  Qn  and  Q22» 

h(x  ,x  )  «  x  +  x 
12  12 

This  is  an  example  of  linear  measurements.  As  we  found  above,  the  RDB  is 

RDB  -  (R-1  +  Qn-1)"1  *  (3-133) 

The  CRVB  is  also  the  minimum  mean  square  error  for  this  case  and  is  given  by 

CRVB  -  (R“l  +  Qn-1  ~  R~2lQ22_1  +  R-1]-1)-1  .  (3-134) 

Thus,  in  this  case  we  always  have  CRVB  >  RDB.  Figure  3-10  shows  the  two 
bounds  versus  R  for  Qn  *  Q22  *  1.0. 

h(x1,x2)  ”  +  *2 


We  computed  the  RDB  above: 
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RDB  -  (UQnlKT1  +  Qll"1)"1  .  (3-135) 

The  CRVB  is  easily  found  to  be 

CRVB  =  ([4QU]R-1  +  Qn-1)-1  .  (3-136) 

Thus,  in  this  case  we  always  have  RDB  >  CRVB.  Figure  3-11  shows  the  two 
bounds  versus  R  for  Qu  =  1.0. 

3. 4. 4. 2  Comparison  at  Low  Signal-to-Noise  Ratios 

We  can  prove  a  general  asymptotic  relationship  between  CRVB  and  RDB  as 
R+°°.  This  relationship  is  similar  to  the  one  we  proved  in  subsection  3. 3. 3. 4. 
Let  x  be  partitioned  into  components  xj  and  X2 ,  and  define  Qn,  Q12*  021*  and 
Q22  in  subsection  3. 4. 2.1.  Recall  that 

Ql  *  Qll  “  Q12  022_1  021  (3-137) 

and 

rl  -  E([h(xi,X2)  -  E(h(xi ,X2) |x2))2)  .  (3-138) 

The  RDB  is  simply 

RDB  -  ((I^Qj-ljR-1  +  Qi-1)-1  .  (3-139) 


The  CRVB  gives  a  lower  bound  on  the  error  covariance  matrix.  This  matrix 
bound  is 


B  -  (Q-1  +  R-1  EdLJ!?  )  r1 
3x  3x 


(3-140) 


The  CRVB  for  the  Xj  estimate  is  the  Bj j  element  of  the  matrix  B.  For  large  R 
we  can  approximate  B  by 
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Figure  3-11.  Comparison  of  Rate  Distortion  and  Cramer-Rao-Van  Trees 
Bounds  for  Vector  State  With  Nonlinear  Measurement 
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B  =*  Q  -  R_1QE( 


3hT  3h 
3x  3x 


)Q  +  0(R-2) 


If  Qj2  *  0,  Q21  *  0,  then  we  have  to  first  order  in  R"*: 

,  3h  „  „ 

CRVB  =  QU  -  R-1Qir2(  [ - ]2)  +  0(R"2) 

3x^ 


Equivalently,  we  have 


CRVB  -  (  [E(  [ - ]2)]r-1  +  QU-1  +  0(R-2))-l 

3xi 


to  compare  with  Eq.  3-139.  From  the  scalar  state  inequality  Eq.  3-76 


,  3h 

Ql-1  *  Ti(x2)  <  E([ - 1 2 | X2 ) 

3xi 

and  consequently, 

Ql-1  Ti  <  E([— ]2) 

3xi 


We  will  assume  that  the  partition  xi ,x?  of  x  has  been  chosen  so  that 
are  independent  (i.e.,  Qi?  »  0,  Q91  ”0). 

Then  we  see  that  Qi  *  Qn  and 

CRVB  <  RDB  +  0(R-2) 


(3-141) 

(3-142) 

(3-143) 

we  have 

(3-144) 

(3-145) 

li  and  x? 

(3-146) 


That  is,  there  is  a  term  0(R“2)  which  converges  to  0  as  fast  as  R“2  when  R+«», 
and  the  CRVB  is  larger  than  RDB  by  at  most  0(R”2).  Note  that  if  h  depends  on 
xi  in  a  nonlinear  way,  then  we  must  have  CRVB  <  RDB  for  sufficiently  large  R. 
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This  is  the  case  in  the  example  h(xi,x2>  ■  xl^+*2  above.  In  the  linear  exam¬ 
ple  h(xj ,X2)*xi+X2 »  one  can  see  from  Eqs.  3-133  and  3-144  that  CRVB  exceeds 
RDB  only  by  a  term  of  order  R”2 

3.4.5  Rate  Distortion  Bound  for  Vector  State  and  Vector  Measurement 

If  an  m  dimensional  vector  measurement  is  taken,  the  results  of  subsection 
3.4.1  change  as  follows.  The  earlier  result  Eq.  3-19  for  the  scalar  state  case 
implies 

E([xi  -  xi(y)]2|x2>  >  Qi  (det(  T] (x2>R_1  +  Im))"*1  (3-147) 

where  Im  is  the  mxm  identify  matrix,  Qj  is  given  as  in  En.  3-99  and 

Ti(x2)  *  E( [h(xi ,X2>  ~  E(h(xi ,X2> |x2)l  [h(xi,X2)  -  E(h(xi ,X2> |x2> ]T |x2> 

(3-148) 

Thus  we  have 

E([X1  -  xi(y)]2)  >Qj  •{E(det(r1(x2)R"1  +  Im))}_1  .  (3-149) 

The  right-hand-side  of  Eq.  3-149  is  computable  in  closed  form  if  each  component 
of  h  belongs  to  the  class  An.  Unfortunately,  this  computation  appears  to  be 
difficult  in  general,  and  further  development  is  required  to  make  the  RDB  use¬ 
ful  for  the  general  vector  state,  vector  measurement  case.  However,  if  the 
measurement  vector  consists  of  N  identically  distributed,  conditionally  inde¬ 
pendent  scalar  measurements,  then  RDB  is  given  by 


RDB  -  (iNrxQi-ijR-1  +  Qi-1)-1 


(3-150) 
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3.5  CONCLUDING  REMARKS 

3.5.1  Summary 

In  this  report  we  have  described  how  to  compute  analytically  rate  distor¬ 
tion  bounds  of  mean  square  error  for  static  nonlinear  estimation  problems  of 
the  form 

y  =  h(x)  +  v  (3-151) 

where  x  and  v  are  Gaussian  distributed.  Specifically,  we  obtained  a  lower 
bound  of 

ECIxi  -  xj (y) ]2)  (3-152) 

where  x^  is  a  scalar  component  of  x.  We  showed  that  the  rate  distortion 
bound  is  asymptotically  tighter  than  the  Cramer-Rao-Van  Trees  bound  in  the 
limit  as  the  noise  covariance  R  becomes  unbounded  (i.e.,  as  signal-to-noise 
ratio  approaches  0).  We  illustrated  the  rate  distortion  bound  and  its  com¬ 
parison  to  the  Cramer-Rao-Van  Trees  bound  using  a  number  of  simple  examples. 

3.5.2  Conclusions 

Based  on  present  results,  the  rate  distortion  bound  offers  a  better 
approximation  of  mean  square  performance  in  the  high  measurement  noise  regime 
than  the  Cramer-Rao-Van  Trees  bound.  Furthermore,  the  rate  distortion  bound 
requires  little,  if  any,  more  computation  than  the  Cramer-Rao-Van  Trees  bound. 
Thus,  the  rate  distortion  bound  appears  to  complement  the  Cramer-Rao-Van  Trees 
method  in  the  nonlinear,  high  noise  regime  where  the  latter  bound  is  known  to 
give  overly  optimistic  approximations  of  the  true  mean  square  error.  However, 
in  order  to  make  the  rate  distortion  theory  useful  for  the  dynamic  nonlinear 
estimation  problems  of  tracking,  we  must  develop  our  current  results  in  two 
significant  ways: 
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1 .  it  is  necessary  to  obtain  a  simple  rate  distortion  bound  in 
the  case  of  a  vector  state  and  a  general  vector  measurement; 

2.  it  is  necessary  to  derive  a  recursively  computable  bound  for 
dynamic  estimation  problems. 

3.5.3  Other  Work 

Zakai  and  Ziv  (11]  first  applied  rate  distortion  theory  to  mean  square 
performance  analysis  of  nonlinear  filtering  problems.  The  results  of  [11] 
were  restricted  to  a  special  class  of  continuous-time  processes.  Galdos  [12] 
extended  these  results  to  general  vector  processes,  both  in  continuous  and 
discrete  time.  In  this  section  we  have  derived  some  preliminary  bounds  on 
individual  component  errors  as  in  Eq.  3-152.  The  results  of  [12]  give  bounds 
on  the  sum  over  all  component  errors 

n 

l  -  xk(y)]2)  .  (3-153) 

k=l 

We  believe  the  approach  here,  based  on  our  earlier  work  [13], [14],  will  yield 
a  more  accurate  estimate  of  mean  square  estimation  error.  However,  until  re¬ 
sults  of  this  section  are  extended,  we  can  make  no  comparisons  with  [11], [12]. 

3. 5. A  Further  Investigation 

Other  directions  for  further  investigation  exist  beside  the  two  necessary 
extensions  noted  in  subsection  3.5.2  above.  One  direction  would  extend  the 
bounds  to  problems  for  which  x  is  an  unknown,  non-random  parameter  (or  a  mix¬ 
ture  of  random  and  non-random  parameters).  In  [1]  we  found  that  a  large  class 
of  tracking  problems  can  be  modeled  by  a  state  process  which  consists  of  an 
unknown  deterministic  component  and  an  unknown,  Gaussian  distributed  random 
component.  Rate  distortion  theory  for  nonstatistical  sources  (e  -  entropy 
methods  [4]  may  allow  us  to  derive  such  results. 
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Another  direction  is  to  study  the  effect  of  architecture  constraints, 
such  as  preprocessing  of  measurements,  on  tracking  estimation  performance. 

We  investigated  this  problem  in  [1]  using  Cramer-Rao-Van  Trees  methods.  A 
rate  distortion  approach,  based  as  it  is  on  info  nation  theory,  would  provide 
a  more  general,  more  accurate  method  of  analyzing  architectural  constraints. 
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SECTION  4 

AMBIGUITY  PERFORMANCE  ANALYSIS 

4.1  INTRODUCTION 

Ambiguity  analysis  ([7], [15]  Chapter  10)  is  an  attempt  to  understand  the 
global  nature  of  a  parameter  estimation  problem.  This  is  in  contrast  to  Cramer- 
Rao  methods  which  provide  a  more  local  analysis  of  estimation  performance. 

The  Cramer-Rao  lower  bound  on  mean  square  estimation  error  will  be  an  accurate 
estimate  of  true  optimal  performance  provided  that  it  is  possible  to  acquire 
or  maintain  an  estimate  near  the  unknown  parameter  at  all.  The  local  problem 
(addressed  by  the  Cramer-Rao  method)  is  to  analyze  accuracy  given  acquisition; 
the  global  problem  (addressed  by  ambiguity  analysis)  is  to  analyze  the  acqui¬ 
sition  performance. 

The  ambiguity  approach  can  be  formulated  as  follows.  Suppose  one  wishes 
to  estimate  an  unknown  parameter  x  given  a  measurement  y.  Let  x  denote  the 
maximum  likelihood  estimator  of  x  which  depends  on  y  and  consider  the  mean 
square  estimation  error 

M(x  '  x)2}  -  l  Ex{(x  -  x)2jxeRk}  Px{xeRk} 
k-0 

where  E^t*}  and  PxW  denote  the  expectation  and  probability  given  that  x  is 

the  true  value  of  the  parameter;*  and  R.,  R,,  ...»  R  are  nt-1  regions  subdivid- 

oi  n 

ing  the  x  parameter  space.  Approximate  Ex{(x  -  x)2|xeR^}  by  e^2  and  approxi¬ 
mate  P^xeR^}  by  p^.  Then  the  mean  square  error  is  approximated  by 
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Ex{(x  -  *)2}  -  l  £k2  *  Pk  . 
k-0 

For  example,  assume  that  the  true  parameter  xeRQ  and  let  eQ2  be  the  Cramer-Rao 
bound  for  the  problem.  If  k  *  0,  choose  a  typical  x^cR^  and  let  q^2  *  (x  -  x^) 
Approximate  Px{xcRk}  by  the  following  hypothesis  testing  problem.  Let  xn(y) 
be  the  x^  that  maximizes  P^fy)  and  lat  Pk  ^xi^n^k}* 

Can  we  show  rigorously  that 

l  ek2  *  Pk  - ►  Ex{(x  -  x)2} 

k-0 

as  the  number  nfl  of  regions  increases  and  the  size  of  the  regions  decreases? 
Can  we  estimate  the  size  of  the  error  for  a  given  n  and  choice  of  regions  R^? 

In  this  section  we  provide  detailed  convergence  analysis  for  a  partic¬ 
ular  sequence  of  approximations  for  the  calculation  of  the  error  variance  in 
a  maximum  likelihood  estimation  problem.  We  restrict  our  attention  here  to  a 
scalar  problem.  While  several  of  the  detailed  calculations  we  perform  do  use 
the  scalar  nature  of  the  problem  to  allow  us  to  write  down  very  explicit  for¬ 
mulae,  the  general  nature  of  the  analysis  can  be  extended  (this  would,  how¬ 
ever,  involve  the  determination  of  several  additional  estimates  to  replace 
the  closed-form  expressions  available  in  the  scalar  case) . 

A. 2  PROBLEM  FORMULATION 

We  consider  the  problem  of  estimating  a  scalar  parameter  x  which  is  known 
to  take  values  on  the  interval  [0,1].  We  have  available  the  scalar  measurement 

y  -  h(x)  +  v  (A-l) 

where  v  is  a  Gaussian  random  variable  with  mean  0  and  variance  1.  We  also 


assume 
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h(0)  >  0 


(4-2a 


0  <  h'(x)  <  M  <  *»  for  all  xe[0,lj 


(4-2b 


and  h(x)  can  be  expanded  In  a  series  around  any  point  ae[0,l] 


h(x)  -  h(a)  +  —  (a)[x  -  a]  +  R(x  -  a) 
dx 


(4-3) 


where  r(x  -  xQ)  ■  o((x  -  xQ)2).  Note  that  Eq.  4-2a  is  a  trivial  assumption 
since  we  can  always  add  a  constant  to  h.  Also  note  that  the  monotonicity 
assumption  (Eq.  4-2b)  simply  avoids  the  possibility  that  h(Xj)  ■  h(x2)  for 
any  x^  x2e[0,l]. 

The  problem  with  which  we  are  concerned  is  the  following.  Suppose  that 
the  true  value  of  x  is  xQ.  We  wish  to  calculate  (or  more  precisely  to  obtain 
a  sequence  of  approximations  to) 


E[(x  -  xj2|x  «*  x  ] 


where  x  is  the  maximum  likelihood  estimate.  That  is,  let 


£(x)  ■  yh(x)  -  h(x)2 

2 


(4-4) 


x  *  arg  max  £(x) 


(4-5) 


We  begin  with  several  preliminary  calculations. 


Computation  of  the  Distribution  Function  for  x 


Let  us  rewrite  Eq.  4-4  using  Eq.  4-1  and  the  fact  that  x  ■  xQ: 


£(x)  ■  h(x.)  h(x)  +  v  h(x)  -  -  h(x)2 


(4-6) 
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Consider  the  derivative  of  Eq.  4-6 


V 

(x)  ■ 

■  [h(xQ)  +  v 

-  h(x)  ]h'( 

x) 

(4- 

■7) 

Note 

that,  thanks  to  Eq. 

4-2b 

r(x)  -  o  if 

and  only 

if 

v  -  h(x)  - 

h(x0) 

(4- 

-8) 

Thus 

,  again  thanks  to  Eq. 

4-2 

we  see  that 

*’(x)  >  0 

for 

all  xe[0,l] 

if  v  >  h(l)  -  h(xQ) 

(4- 

-9a) 

i'(x)  <  0 

for 

all  xe[0 ,1 ] 

if  v  <  h(0)  -  h(xQ) 

(4- 

-9b) 

Thus 

Prob(x  = 

■  0|x 

=  xQ )  "  Proh 

.(v  <  h(0) 

-  h<x0>) 

(4- 

-10a) 

Prob(x  - 

■  1  j  X 

“  xQ)  -  Prob 

»(v  >  h(l) 

-  h(xQ)) 

(4- 

-10b) 

If  ve(h(0)  -  h(xQ),  h(l)  -  h(xQ)),  x  will  be  the  value  of  x  for  which  l' (x)  *  0, 
i.e.,  the  value  for  which  Eq.  4-8  is  satisfied.  Thus, 

Prob(0  <  x  <  a|x  »  xQ)  -  Prob(h(0)  -  h(xQ)  <  v  <  h(  a)  -  h(xQ)) 

/•h(  a)-h(x  )  (4-11) 

■  /  N(v;  0,1)  dv 

h(0)-h(xg) 

where 

N(v;  0,1)  -  — l—  e"v2/2  (4-12) 

V2i 


Thus  the  probability  density  function  for  x  on  [0,1]  is 
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P„(a)  »  —  Prob(0  <  x  <  a|x  -  xft) 
x  da  u 

=  N(h(a)  -  h(xQ);  0,l)h'(a) 


(4-13) 


Note  that  in  this  case  we  obtain  a  formula  for  the  error  variance 
e[(x  -  xQ)2|x  *>  xQ]  =  xQ2  Prob(v  <  h(0)  -  h(xQ)) 

+  (l  -  x Q ) 2  Prob(v  >  h(l)  -  h(xQ)) 

r  1  h*(a)  (  !  ) 

+J  (xQ  -  a)2  ~^ZT  exp  j  -  -  [h(  a)  -  h(xQ)  )2  j  da 


In  the  sequel  we  develop  a  sequence  of  approximations  to  this  quantity  moti¬ 
vated  by  our  desire  to  develop  methods  that  can  be  applied  to  more  complex 
problems. 


The  Cramer-Rao  Bound  and  the  Ambiguity  Function 
For  this  problem 


p(y|x)  =  - -  exp 

yfz i? 


1 

2 


(y  -  h(x)  )2 


(4-14) 


and  it  is  a  straightforward  computation  to  verify  that  the  Cramer-Rao  Bound  is 


CRB  (xQ) 


H32  in  p(y|x  ) 


(4-15) 


The  ambiguity  function  in  this  problem  is 


A(x  ,x^)  “  h(x1)h(x2) 


(4-16) 


4-5 
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Note  that 


e|^£<x)|x0-  x]  -  A(x,xq)  -  i  A(x,x) 
Cov  [i(x1),£(x2)  ]  -  A(xlfx2) 


(4-17) 

(4-18) 


and 


il 

3x2 


|^A(x,x0)  -  -  A(x,x)^j 


1-1 


CRB  (xQ) 


(4-19) 


x  *  x. 


Thus  the  Cramer-Rao  bound  Is  seen  to  depend  explicitly  on  the  curvature  of 
E[£(x)|x  ■  xQj  at  the  location  of  Its  peak,  i.e.,  at  x  -  xQ. 


4.3  CONVERGENCE  ANALYSIS 


We  now  construct  a  sequence  of  approximations  indexed  by  the  Integer  N. 
Essentially  what  we  will  do  is  to  divide  the  interval  [0,1]  into  subintervals, 
most  of  which  will  be  of  length  1/N.  There  will,  however,  be  one  interval 
centered  at  the  true  value  xQ  that  will  be  larger.  Specifically,  let 


xo  + 


1 

vr 


n  lo,i] 


(4-20) 


Assuming  that  N  is  large  enough  so  that  xn  -  1/  >  0,  define 


Il,-1N  -  (0) 


(4-21a) 


L(N) 


(4-21 b) 


where  a  A  b  Indicates  minimum  of  a  and  b  and  where 
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L(N) 


-yJW -  1 


(Here  ph  is  the  smallest  integer  greater  than  or  equal  z) . 
L(N)  -  0(N) .  Similarly  assuming  that  xQ  +  1/Vn”<  1»  define 


(4-21c) 


Note  that 


Ir,-in  -  {1} 


(4-22a) 


IR»i 


N  . 


,  i  -  0,  1, 


where  a  V  b  indicates  maximum  of  a  and  b  and  where 


R(N)  = 


R(N) 

(4-22b) 


(4-22c) 


Again  R(N)  ”  0(N).  What  we  have  done  is  to  partition  [0,1]  into  disjoint 
sets.  There  is  one,  larger  central  set  ICN,  and  the  two  endpoints  Il,-1N  and 
The  remaining  sets  to  the  left  of  ICN  are  of  the  form  (i/N,  i+l/N] 
except  for  the  one  bordering  on  ICN  which  is  clipped  off  so  that  it  doesn't 

overlap.  Similarly  the  sets  to  the  right  of  ICN  are  of  the  form  [i/N,  i+l/N) 

except  for  the  one  bordering  on  Ic^  which  is  clipped  off  so  that  it  doesn't 

overlap. 

Since  these  sets  don't  overlap,  we  have  the  following  equality 

L(N) 

E[(x  -  x0)2la  "  x0]  ”  l  EUX  “  X0)ZIX  “  xo»  xelLiNMxelLiNlx  *  x0) 

i— 1 
R(N) 

+  l  E[(x  -  X0)2IX  "  x0»  X£lRiN]Pr(xelRilx  M  x0) 
i— 1 

+  E((x  "  X0)2IX  “  xo»  xeIcN]Pr(jJeIcN|x  -  x  ) 

0  (4-23) 


4-7 
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I  I 

i  - 

i 


1 1 


i 


i 


Our  approximate  method  for  evaluating  the  left-hand  side  of  Eq.  4-23  is  based 
on  obtaining  approximations  for  each  of  the  terms  on  the  right-hand  side.  To 
do  this,  we  proceed  by  defining  the  following  discrete  set  of  points 

-  0 

5^N  «  center  point  of  Il,in  >  i  =  0,  L(N) 

xo  (4-24) 

Y~lN  *  1 

Yi1*  ■  center  point  of  Ir^N  ,  i  =  0,  ...»  R(N) 


Note  that 


60N  - 

61+1N  -  6iN 

5i+1N  -  5iN 
*0  ~  $L(N)N 

yr(n)n  -  *0 

YiN  "  Yi+1N 
YiN  -  Y1+1N 


1 

2N 

1 

N 


•(:) 


l 

N 

1_ 

2N 


,  i  “  0,  . L(N)  —  2 

,  i  =  L(N)  -  1 


,  i  -  R(N)  -  1 

,  i  -  0,  R(N)  -  ° 


(4-25) 


4-8 


I 
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Consider  Che  hypothesis  testing  problem  in  which  we  assume  that  we  know 
that  x  takes  on  one  of  the  finite  set  of  values  in  Eq.  4-24,  and  suppose  that 
we  use  Che  maximum  likelihood  criterion  for  choosing  our  estimate  from  this 
finite  set.  Let 

PLiN(xQ)  -  Prob(choose  6iN|x  -  xQ)  ,  i  -  -1,  L(N)  (4-26a) 

PCN(X0)  “  Prob( choose  xQ|x  -  xQ)  (4-26b) 

PRiN(x0)  "  Prob( choose  YiN|x  -  xQ)  ,  i  -  -1,  ....  R(N)  (4-26c) 

Then  our  approximation  to  Eq.  4-23  is 

T  «  "1  L<N> 

E[^(x  -  x0)2|x  ,  XqJ  PliN(x0) 

R(N) 

+  JL/V  "  X0)2  PRiN(X0}  (4‘27) 

+  CRB  (x0)p(N)  PCN(X0)  &  v(xQ) 

where 

P(N)  =  E  v2 1  J  v|  <  — —  —  (xn)  (4-28) 

LI  v*r  dx  °  J 

We  now  proceed  to  estimate  the  errors  in  the  various  terms  and  to  show  that 
Eq.  4-27  converges  to  e[(x  -  xQj2|x  =  x^j  as  N 


♦  <*>. 
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The  Terms  b[(x  -  xQ)2|x  -  xQ>  xeIL>1N]  and  e[(x  -  xQ)2[x  -  xQ>  xeIR>iN] 

Note  first  that 

E[(x  ~  X0)2IX  "  x0.  xeIL,-lN]  "  ($-lN  "  x0)2  (4-29a) 

E[(x  "  X0)2IX  "  x0»  *eIR,-lN]  "  (y-lN  ~  x0)2  (4-29b) 


Furthermore  for  i  >  0 


E[(x  ~  xn)2lx  *  xo»  xeILiN] 


(4-30) 


f  (“  “  x0)2  P4“lx  *  x0»  XElLiN)  dt 


I  N 
Li 


where 


Pi(°>lx  -  x0,  xeILiN) 


Px(“l 


x  “  xn) 


J‘  P-(a|x  -  x  ) 

*'lLiN  xV  0 


(4-31) 


and  p»(ot|x  ■  x.)  is  given  in  Eq.  4-13.  We  now  see  that  the  integral  in  Eq. 
x  u 

4-30  is  sufficiently  smooth  over  the  entire  interval,  so  that  we  may  approxi¬ 
mate  it  as 


(«iN  -  xQ)2  p«(6iN|x  -  xQ,  x£Ilin)  *  Length (lLiN)  (4-32) 


The  error  introduced  in  this  approximation  is  0 
Fig.  4-1. 


To  see  this,  examine 
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Figure  4-1. 


Approximation  Error 


The  integral  to  be  evaluated  is  the  area  under  the  solid  curve ,  while  Eq. 
4-32  equals  the  area  under  the  dotted  curve.  Let 


I  ^  f"”  * 

K  -  sup  sup  —  (a  -  x  )2  P«(a|x  =  x  ,  xeILiN) 
i  a  |  do  L  u  x 

It  is  easy  to  see  that  K  <  »  and  that  the  magnitude  of  the  difference  between 
the  areas  under  the  solid  and  dashed  curves  in  Fig.  4-1  is  bounded  above  by 


K 

2 


Length (lL1N) 


I 


2 


which,  from  Eq. 
argument 


4-21  is  0 


Finally  we  note  that  by  the  same  type  of 


1 


mJ  p-Nx  - 

*Li 

-  P -( <5iN|x  -  x  , 
x  0 


0,  xeILiN)  da 

xeILiN)  *  Length (lLiN)  +  0^  —  j 


(4-33) 
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Combining  Eqs.  4-30  through  4-33  we  see  that 


Ef(*  -  x0)2|x  =  xQ,  xeILiN]  =  (6iN 
i  -  0,  ....  L(N) 


In  an  analogous  fashion  we  can  show  that 


E[(*  -  xQ)2|x  -  x0,  xeIR1N]  =  (Y1N 
1=0,  ...»  R(N) 


The  Term  e[(x  -  xj2|x  *  x„,  xeIcN] 


Substituting  Eq.  4-3  into  Eq.  4-6  we  obtain 


1  ,  1  fdh  ,  2  dh 

l(x)  =  -  h(xfl)2  -  -  (xQ)  J  [x  -  xfl]  -  -  (xQ)[x  -  xq]r(x  -  xQ) 


-5[ 


~|2  dh 

R(x  -  xq)  +  v  h(x  )  +  v  -*  [x  -  x  ]  +  v  R(x  -  x  ) 


(4-36) 


Assuming  that  xeIcN  we  have  that 


d*(x) 


fdh  ,  ,  "y2  ,  dh  , 

~  Is  W  J  [x  "  XoJ  +  V  Tx  W  +  A1(X)  +  VA2(X) 


(4-37) 


where 


=  •  £  (X0)[R(^  -  *„)  +  (*  -  x0)  ^  R(x  -  x0)] 
-  R(x  -  x  )  —  R(x  -  x J 


=  —  Rfx  -  x  ) 


(4-38a) 


(4-38b) 
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and 


sup 

xel  N 


c 


Aj(x)| 


sup 
xel  N 


c 


A2(x)} 


(4-39a) 


(4-39b) 


From  Eq.  4-37  we  see  Chat 


x  -  x„ 


dh  (  rr1  rdh  ,  n~2r  -  -n 
z  (*0)J  v+  [z  WJ  [_Ai(x)  +  vVx) j 


(4-40) 


Also  we  are  assuming 


Combining  Eqs.  4-39  and  Eq.  4-40  we  can  deduce  that  the  Implied  constraint 
on  v  is 


Thus 


(4-41) 


E  v2  |x  *  x.  ,  xel, 


,»]-«[^||v|  <-^=- 5  («,)]+ 
°(^rJ 


(4-42) 


P(N)  + 


The 


VN5/2  J 


comes  from  the 


•(=) 


terra  in  Eq.  4-41  which  implies  that  the 
actual  limits  on  v  can  differ  by  a  term  of  order  1/N.  Thus  the  probability 
mass  in  the  interval  between 


4-1 
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1  dh 

yJT  dx 


and  the  right-hand  side  of  Eq.  4-41  is  0 


and 


1  dh  m"l2  f  1  dh 

YiT  dx  0  'N  /  j  yjr  dx 


Finally,  using  Eqs.  4-39,  4-40,  and  4-42  we  obtain 


xel 


P(N) 


“dh 

_dx 


(XJ 


iT  +  of-i-) 

J  VN3/2>I 


(4-43) 


(4-44) 


The  Probabilities  in  Eq.  4-23 

Under  the  assumption  that  x  is  one  of  the  points  in  Eq.  4-23,  the  maximum 
likelihood  decision  rule  is  to  choose  x  corresponding  to  the  largest  among  the 
values  A(x)  evaluated  at  these  points.  Note  next  that  under  the  assumption 
that  x0=*  x 

*(di+lN)  _  *(diN)  “  h(x0)j^h(,si+lN)  “  h(6iN)^J 

-  |  £h(6i+1N)2  -  h(6iN)J2  (4-45) 

+  v£h(5i+1N)  -  h(6iN)^| 

Thus  using  Eq.  4-2b  we  see  that  the  sign  of  -  l(di^)  is  the  same  as 

the  sign  of 


4-14 
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hK)  -  \  [  h(si+lN)  +  h(6iN)  1  +  v  ^  ~WLiN  +  v  (4-46a) 

Similarly  the  signs  of  i(xQ)  -  *(«l(N)N)»  *(^R(N)N)  "  *(x0)>  and  *(TLN) 

-  *(n+iN)  are  the  same  as  the  signs  of 


h(x0) 

1 

2 

h(xQ)  +  h(«L(N)N) 

+  v  A  ~*1clN  +  v 

(4-46b) 

o 

X 

Jd 

1 

2 

[^h(YR(N)N)  +  h(x0) 

+  v  ^  -pc2N  +  v 

(4-46c) 

h(x0) 

1 

2 

h(fiN)  +  h(TL+lN) 

+  v  ^  “MRi^  +  v 

(4-46d) 

Again  using  Eq.  4-2b)  we  have  that 

PL,-1N  <  WL,0N  <  •••  <  PL,L(N)-1N  <  PclN  <  >jc2N  <  PR,R(N)-1N  <  •••  <  Pr,-1N 

(4-46e) 

and  from  this  we  can  deduce  that  the  quantities  in  Eq.  4-46  are  as  follows: 


pl,-in(x0) 

-  Prob( v  <  Pl,-1N) 

(4-47 a) 

PL,iN(xo) 

-  Prob(yL>i_1N  <  v  < 

i  -  0,  1 . L(N)  -  1 

(4-47b) 

PL,L(N)N(x0) 

-  Prob(pL>L(N)_1N  <  v  <  pclN) 

(4-47c) 

PCN(X0) 

*  Prob(uciN  <  v  <  pc2N) 

(4-47d) 

pr,r(n)n 

-  Prob(pc2N  <  v<  Pr,r(n)-1N) 

(4-47e) 
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PR,i 


Prob(wRflN  <  v  <  HR,i-iN) 
i  =  0,  R(N)  -  1 


PR,-1  “  “  Prob( v  >  Pr,-in) 


(4-47f ) 


(4-47g) 


We  now  compare  these  terms  to  the  terms  to  which  they  correspond  in  Eq. 
4-23.  First,  note  that 

PL,-1N(*o)  =  Prob( choose  0|x  -  xQ)  -  Prob  ^  v  <  i  +  h(0>^J  “  h^X0^) 


(4-48a) 


while 


Prob(xelL>_iN|x  *  x  )  *  Prob(v  <  h(0)  -  h(x  ))  (4-48b) 


Given  Eq.  4-2b  we  have 


Similarly 


Prob(xelLf-1N|x  -  x0)  -  p^N  +  0^  ) 


Prob(xeIR>-1N|x  -  x  )  =  PR,-!1*  + 


(4-49) 


(4-50) 


Next  note  that  for  i  -  0,  1,  ...,  L(N)  -  2 


-  J  N(v;  0,1)  dv  “  J"  N(v;  0,l)  dv 


PL,i-l 


N 


) 


(4-51) 
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while  from  Eq.  4-11  for  i  -  0,  1,  L(N)  -  2 


;(tK> 

Prob(xelL>iN |x  =  xQ)  =  /  N(  v;  0,1)  dv 

h(;)'h(x») 


(4-52) 


Using  Eq.  4-3  (with  a  =  i/N)  we  have  that 


(4-53a) 


if  (l  +  3/2\  /i  +  l/2'i~]  f  i  +  1  /I  ^ 

i  LhHr~J  +  hrr-JJ  J  +  0bJ  <4'53b) 


Prob(xeIL>iN|x  -  xQ)  -  PL,iN(xQ)  +  °(^)  ,  i  »  0,  1 . L(N)  -  2 


(4-54) 


Similarly 


Prob(iclRtiN|x  -  xQ)  -  PR,iN(x0)  +  O^J  ,  i  -  0,  1 . R(N)  -  2 


(4-55) 


Now  for  i  -  L(N)  -  1 


\  [h(6L(N)N)+h(--^"—  )]]-h(x0) 
pl,l(n)-in(*0)  -  /  N(v;  0,1)  dv 


(4-56a) 


i 
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while 


/ L(N)  'S  . 

h(.T'J'h(x«) 

Prob(xeIL>L({j)_1N|x  =  *Q)  =  /  N(  v;  0,1)  dv 

('L(N)-n  . 

hl~J‘h(l,o) 


(4-56b) 


and  furthermore 


{l<«n  •  1  [  ir  +  ( x«  -  vH] 


(4-57) 


Using  Eq.  4-3  we  have 


so  that 


Prob(xelL ,L(N)-1N I x  -  x0)  *  PL,L(N)-1N(*0)  +  0 [ Jj )  (*”59) 


and  similarly 


prob(ieIR>R(N).1N|x  -  xQ)  -  PR,R(N)-1  (x„)  +  °(")  (4-60) 


Next  we  have  that 


”  ^h(x0)+h(SL(N)N)^]-h(xo^ 

pl,l(n)n(x0)  -/  N(v;  °»i>  dv 

2  [h(6UN)N)-H  h  ^  ]  J-h(xQ) 


(4-61) 


4-18 


1 
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and 


»k-  — 

>•  0  vr 

Prob(xeIL>L(N)N)|x  =  xQ)  =  /  N(v;  0,1)  dv 


A  similar  argument  to  the  preceding  ones  yields 


Prob(ieIL>L(N)N|x  -  xQ) 


pl,l(n)n 


and  similarly 


Prob(xeIR>R(N)N|x 


x0)  =  PR,R(N)N 


Finally 


PcN(*0) 


"  j^h(x0)+h(^(N)N)^]~h(> 


J  N(v;  0,1)  dv 

\  ^h(x0)+h(6UN)N)^|-l'(x0 


r 
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Using  Eq.  4-57,  the  corresponding  expression  for  YR(N)N  flnd  Eq.  4-3  we  find 
that 


i  [h(x»)  +  >]  ■  h(x°  •  f  °(^f) 


i  [htx°>  +  h(YM»>N  >]  •  h(x0  +  ^=0  +  °(^=-) 


(4-67) 


(4-68) 


and  thus 


Prob(xeIcN|x  -  xQ)  -  pcN(xq)  +  0  (~r) 


(4-69) 


Combining  the  estimates  4-35,  4-44,  4-49,  4-50,  4-54,  4-55,  4-59,  4-60, 
4-63,  4-64,  4-69,  and  the  facts  that  P(N)  =*  of -  |  and 

^  Vn~  ' 


{- 


E  (x  “  X0)2IX  *  xq»x£IL(N) 


£(i  -  x0)2|x  „  x0,ieIR(N)N^j  - 


(4-70a) 


(4-70b) 


we  conclude  that 


£(x  -  x0)2|Xq  -  cTJ  «  v(xQ)  +  0  j 


(4-71) 


where  v(x0)  is  the  approximation  in  Eq.  4-27. 


4-20 


9  J 
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4.4  CONCLUSION 

In  this  section  we  have  studied  the  ambiguity  analysis  method  of  approx¬ 


imating  mean  square  estimation  error  for  a  simple  nonlinear  parameter  estima¬ 
tion  problem  given  by  Eq.  4-1.  We  saw  that  the  method  outlined  in  subsection 
4.1  converges  to  the  true  mean  square  error  of  the  maximum  likelihood  estima¬ 
tor  and  that  the  error  is  inversely  proportional  to  the  number  of  regions  used 
to  subdivide  the  parameter  space.  Note  that  the  central  region  RQ  for  which 
the  Cramer-Rao  bound  was  used  to  estimate  E* {(xQ-x )2|xeRQ }  was  proportional 
to  (vr-i  )  in  size.  The  remaining  regions  were  proportional  to  r1  in  size. 


For  large  N  this  means  that  N  is  of  the  order  of  the  number  n  of  regions 


subdividing  the  parameter  space.  Thus,  the  approximation  error 


also  0 


(:)■ 


•(;) 


Further  work  is  required  to  determine  how  the  magnitude  of  the  measure¬ 
ment  noise  (or  equivalently,  the  signal-to-noise  ratio)  enters  into  the  approx¬ 
imation  error.  This  result  will  clarify  how  large  the  number  of  regions  needs 
to  be  for  a  given  signal-to-noise  ratio.  The  convergence  analysis  also  needs 
to  be  extended  to  the  general  case  of  vector  states  and  measurements  and  to 
the  case  where  a  measurement  process  is  observed.  The  order  of  the  approxima¬ 
tion  error  is  expected  to  remain  the  same  in  these  generalizations  but  a  more 
precise  idea  of  the  size  of  this  error  would  help  us  understand  the  computa¬ 
tional  feasibility  of  applying  the  ambiguity  analysis  method  to  analyze  the 
performance  of  complex  estimation  problems. 
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SECTION  5 

CONCLUDING  REMARKS 


5.1  GENERAL  SUMMARY 

The  research  described  in  this  report  has  investigated  methods  for  pre¬ 
dicting  performance  in  passive  tracking  problems.  Our  objective  has  been  to 
develop  performance  prediction  methods  that  are  computationally  efficient, 
applicable  to  realistic  passive  tracking  models,  and  accurate.  In  our  pre¬ 
vious  work  [1]  we  developed  a  Cramer-Rao  method  to  obtain  a  method  that  was 
computationally  efficient  and  applicable  to  a  large  class  of  mathematical 
models.  In  Section  2  of  this  report  we  have  shown  that  this  method  is  easy 
to  apply  to  more  realistic  models  than  the  ones  used  in  [1],  Specifically, 
we  have  used  the  method  to  study  the  effect  of  uncertain,  unstable  source 
frequency  and  the  effect  of  the  presence  of  a  broadband  source  component  on 
tracking  accuracy. 

In  some  nonlinear  estimation  problems  of  low  signal-to-noise  ratio, 
Cramer-Rao  methods  may  predict  performance  much  better  than  the  optimal  pro¬ 
cessing  algorithm  can  actually  achieve.  This  disadvantage  of  Cramer-Rao  meth¬ 
ods  motivated  us  to  investigate  performance  prediction  methods  which  would 
be  more  accurate  when  the  signal-to-noise  ratio  was  low,  but  which  are  still 
efficiently  computable  for  a  large  class  of  realistic  models.  Sections  3  and 
4  focused  on  this  problem. 
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Section  3  investigated  an  analytical  ( i.e not  requiring  simulation) 
method  based  on  rate  distortion  theory  [4).  This  method  shows  great  promise 
because  it  is  efficient  to  compute  for  a  large  class  of  nonlinear  problems 
and  it  is  better  than  the  Cramer-Rao  method  when  signal-to-noise  ratio  is 
low.  However,  the  method  requires  further  development  to  make  it  applicable 
to  realistic  dynamic  problems. 

Section  5  investigated  a  numerical  method  of  performance  prediction 
often  described  as  ambiguity  analysis  17],[15].  This  method  is  essentially 
based  on  numerical  computations  rather  than  on  analytical  formulas.  The 
method  can  give  an  accurate  performance  prediction  provided  sufficient  com¬ 
putational  resources  are  available.  Our  investigation  studied  the  relation¬ 
ship  between  prediction  accuracy  and  computational  complexity  for  this  method. 
Further  work  remains  to  determine  the  precise  effect  of  signal-to-noise  ratio 
on  the  relationship  between  prediction  accuracy  and  computational  complexity. 

5.2  CRAMER-RAO  PERFORMANCE  ANALYSIS  OF  FREQUENCY  INSTABILITY  AND 

BROADBAND  SIGNALS 

The  method  developed  in  [1]  was  used  to  study  two  effects:  the  effect 
of  an  initially  uncertain  and  randomly  unstable  source  frequency,  and  the 
effect  of  a  broadband  component  of  the  source  frequency.  The  random  fre¬ 
quency  was  parameterized  by  two  variables,  initial  root-mean-square  uncer¬ 
tainty  and  rate  of  variation.  The  rate  of  variation  (studied  for  .01  to  1.0 
Hz/min)  appeared  to  have  negligible  effect  on  both  position  and  velocity 
tracking  error.  Initial  uncertainty  had  little  effect  on  position  tracking 
error,  but  it  did  have  a  significant  effect  on  velocity  tracking  error.  When 
the  initial  source  frequency  uncertainty  reaches  1  Hz,  the  initial  velocity 
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Cracking  error  is  not  substantially  reduced  until  the  source  passes  through 
the  sonobuoy  field.  This  indicates  that  initial  uncertainty  concerning  source 
frequency  can  make  the  velocity  tracking  performance  sensitive  to  source¬ 
sensor  geometry  (i.e.,  good  velocity  tracking  will  depend  more  crucially  on 
good  geometry). 

The  broadband  source  component  was  modeled  as  a  simple  stationary  first- 
order  Markov  process.  The  bandwidth  of  this  process  was  fixed  at  200  Hz 
and  the  ratio  of  broadband  to  narrowband  power  was  varied  from  10-5  to  101*. 
Increasing  this  ratio  increases  the  total  source  power,  and  the  position  and 
velocity  tracking  error  decrease  as  a  consequence.  The  decrease  in  tracking 
error  is  comparable  to  the  decrease  in  error  with  increased  narrowband  source 
signal  power  studied  in  [1].  This  indicates  that  broadband  source  energy  is 
comparable  in  value  to  narrowband  source  energy,  and  that  the  broadband  com¬ 
ponent  of  the  source  signal  can  be  profitably  exploited  by  an  acoustic  signal 
processing  system. 

Many  other  realistic  models  can  be  analyzed  using  the  methods  of  Section 
2.  However,  before  analyzing  the  performance  of  other  realistic  models  of 
passive  acoustic  tracking  problems,  we  need  to  determine  the  degree  of  opti¬ 
mism  inherent  in  the  performance  prediction  of  Section  2.  Sections  3  and  A 
present  one  approach  to  doing  this,  namely  by  trying  to  develop  more  accurate 
performance  predictions  with  which  to  compare  the  methods  of  Section  2.  It 
is  also  desirable  to  compare  these  performance  predictions  to  the  performance 
of  actual  algorithms.  The  methods  of  Section  2  suggest  a  processing  algorithm 
architecture  that  might  realize  the  performance  prediction  in  some  cases  (see 
[1]  for  discussion). 
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5.3  RATE  DISTORTION  PERFORMANCE  ANALYSIS 

In  Section  3  we  showed  how  to  compute  analytically  rate  distortion  bounds 
of  mean  square  error  for  static  nonlinear  estimation  problems  with  a  Gaussian 
distributed  state  vector  and  additive  Gaussian  noise.  Specifically,  we  ob¬ 
tained  a  lower  bound  of  the  mean  square  estimation  error  for  any  specified 
component  of  the  state  vector.  We  showed  that  the  rate  distortion  bound  is 
asymptotically  tighter  than  the  Cramer-Rao  bound  in  the  limit  of  low  signal- 
to-noise  ratio. 

Based  on  these  results  we  conclude  that  rate  distortion  offers  a  better 
approximation  of  mean  square  performance  in  the  low  signal-to-noise  regime 
than  the  Cramer-Rao  bound.  Furthermore,  our  rate  distortion  bound  requires 
little,  if  any,  more  computation  than  the  Cramer-Rao  bound  for  the  special 
class  of  estimation  problems  of  interest.  Thus,  the  rate  distortion  bound 
appears  to  complement  the  Cramer-Rao  method  in  the  nonlinear,  low  signal- 
to-noise  ratio  cases  where  the  latter  bound  is  believed  to  give  overly  opti¬ 
mistic  performance  predictions. 

In  order  to  make  the  rate  distortion  bound  useful  for  dynamic  nonlinear 
estimation  problems  of  tracking,  we  must  develop  our  current  results  to  sim¬ 
plify  the  computations  for  large  dimensional  state  and  measurement  vectors  and 
to  obtain  recursively  computable  formulas  for  dynamic  estimation  problems. 

5.4  AMBIGUITY  PERFORMANCE  ANALYSIS 

In  Section  4  we  analyzed  the  mean  square  parameter  estimation  error  of 
the  maximum  likelihood  method  using  the  ambiguity  analysis  method.  We  showed 
that  this  numerical  method  converges  to  the  exact  mean  square  error  as  the 
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number  of  discretization  regions  increases,  and  the  approximation  error  is 
inversely  proportional  to  the  number  of  regions  used  to  subdivide  the  param¬ 
eter  space. 

Further  work  is  required  to  determine  how  the  signal-to-noise  ratio  quan¬ 
titatively  affects  the  approximation  error.  This  result  would  clarify  how 
large  the  number  of  discretization  regions  need  to  be  for  a  given  signal-to- 
noise  ratio.  The  convergence  analysis  also  needs  to  be  extended  to  the  general 
case  of  vector  states  and  measurements  and  to  the  case  where  a  measurement  pro¬ 
cess  is  observed.  The  order  of  the  approximation  error  is  expected  to  remain 
the  same  in  these  generalizations,  but  a  more  precise  idea  of  the  size  of  this 
error  would  help  us  understand  the  computational  feasibility  of  applying  the 
ambiguity  analysis  method  to  estimation  problems. 
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